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Executive  Summary 


This  document  is  the  compilation  of  the  analysis  and  results  of  the  research  contract  w7714- 
115195/A  entitled  Research  &  development  for  Long  Term  Evolution  (LTE)  wireless  location 
based  on  synthetic  array.  The  report  is  a  compilation  of  the  research  activities  and  outcomes 
for  the  period  of  January  2012  to  March  31  2014. 

The  report  consists  of  two  parts.  This  document  is  part  A  and  part  B  is  a  separate  document. 
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Acronyms 


AN 

Access  node  (eNodeB) 

AP 

Anchor  point 

ANan 

Access  node  that  is  an  anchor  node 

ANfp 

Access  node  that  is  feature  point 

AWGN 

Additive  White  Gaussian  Noise 

AOA 

Angle  of  Arrival 

BCRLB 

Bayesian  Cramer-Rao  Lower  Bound 

BF 

Bayesian  Filter 

BFI 

Bayesian  Fisher  Information 

BFIM 

Bayesian  Fisher  Information  Matrix 

CDMA 

Code  Division  Multiple  Access 

CML 

Concurrent  Mapping  and  Localization 

CRLB 

Cramer-Rao  Lower  Bound 

CV 

Computer  Vision 

EKF 

Extended  Kalman  Filter 

FastSLAM 

Factorized  Solution  for  Simultaneous  Localization  and  Mapping 

FI 

Fisher  Information 

FIM 

Fisher  Information  Matrix 

FP 

Feature  Point 

GNSS 

Global  Navigation  Satellite  Systems 

GPS 

Global  Positioning  System 

i.i.d 

Independent  and  Identically  Distributed 

IMM 

Interacting  Multiple  Model 

IMU 

Inertial  Measurement  Unit 
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IEEE 

Institute  of  Electrical  and  Electronics  Engineers 

KF 

Kalman  Filter 

KL 

Kullback-Leibler 

LAMBDA 

Least-squares  Ambiguity  Decorrelation  Adjustment 

LJG 

Linear  and  Jointly  Gaussian 

LML 

Local  Maximum  Likelihood 

LOS 

Line  of  Sight 

LTE 

Long  Term  Evolution 

MAP 

Maximum  A  Posterior 

MC 

Monte  Carlo 

M-CRLB 

Modified  Cramer-Rao  Lower  Bound 

MEMS 

Micro  Electro-Mechanical  Systems 

Ml 

Mutual  Information 

MLE 

Maximum  Likelihood  Estimator 

MM 

Multiple  Model 

MMSE 

Minimum  Mean  Square  Error 

MN 

Mobile  Node 

NEES 

Normalized  Estimation  Error  Squared 

OCXO 

Oven-Controlled  Crystal  Oscillator 

OWLS 

Opportunistic  Wireless  Localization  System 

PDF 

Probability  Density  Function 

PDOA 

Phase  Difference  of  Arrival 

PEB 

Position  Error  Bound 

PF 

Particle  Filter 

PRN 

Pseudo-Random  Noise 
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POA 

Phase  of  Arrival 

RB 

Rao-Blackwellized 

RBPF 

Rao-Blackwellized  Particle  Filter 

RFID 

Radio  Frequency  Identification 

RIM 

Research  in  Motion 

RTOF 

Round-trip  Time  of  Flight 

SLAM 

Simultaneous  Localization  and  Mapping 

SNR 

Signal-to-Noise  Ratio 

SS 

Signal  Strength 

TCXO 

Temperature  Compensated  Crystal  Oscillator 

TDOA 

Time  Difference  of  Arrival 

TOA 

Time  of  Arrival 

UWB 

Ultra  Wideband 

WiFi 

Wireless  Fidelity 

WLAN 

Wireless  Local  Access  Network 

Variables  used 


variable 

description 

beliX0:t) 

Posterior  pdf  of  the  particles  over  the  complete  time  interval  from  0  to  tTs 

h, 

range  offset 

Bs 

signal  bandwidth 

c 

propagation  velocity 

eixm 

lxm  row  vector  whose  elements  are  zero,  except  the  i  -th  element  which  is 

equal  to  one. 
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/ 

carrier  frequency 

fb 

clock  frequency  bias 

f 

clock  frequency  drift 

m 

clock  random  frequency  error 

I(X,Y) 

mutual  information  between  Y  and  X 

J 

Fisher  information  matrix  in  general 

J tot 

total  Bayesian  Fisher  information  matrix 

Jz 

measurement  information  matrix 

JP 

apriori  information  matrix 

kan 

vector  of  direction 

m 

stacked  vector  of  AN  locations 

m, 

state  vector  describing  the  location  of  the  /th  ANs 

m 

stacked  vector  of  APs  location 

m 

stacked  vector  of  FPs  location 

2D  location  variables  of  an  AN 

NAN 

number  of  ANs 

nap 

number  of  APs 

nfp 

number  of  FPs 

Nd 

dimension  of  dynamic  variable  in  state  vector 

Ns 

dimension  of  stationary  variable  in  state  vector 

NP 

number  of  particles 

Pt 

MN  location  vector 

Pi, 

history  of  MN  locations 

Qi 

state  vector  at  time  step  t 
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dynamic  portion  of  state  vector  at  time  step  t 

<L 

stationary  portion  of  state  vector 

a 

covariance  matrix  of  dynamic  variables  update  process 

motion  process  covariance  matrix 

r(t) 

geometric  range  between  the  AN  and  MN 

s(t) 

bandpass  signal 

s(t) 

lowpass  signal 

u, 

control  vector 

ui* 

history  of  update  control  inputs 

V, 

motion  update  process  noise 

Wut 

measurement  noise 

z, 

measurement  vector  received  at  time  step  t 

Zu, 

observation  vector  from  /'- th  AN  at  the  time  t 

ZW 

history  of  observation  from  /-th  AN 

SiJ 

uk+ 1 

(k  + 1)  x  (k  + 1)  dimensional  matrix  whose  elements  are  all  zero  except  at  the  i  - 

th  row  and  the  j  -th  column  which  is  one 

fit MN 

MN  clock  bias 

St  AN 

AN  clock  bias 

^ MN 

MN  clock  rate  drift 

AN 

AN  clock  rate  drift 

&<p(t) 

carrier  phase  variation 

a 

V N p  ,  second  order  derivative  operator 

carrier  phase  measurement  noise 

V 

normalizing  constant 
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A 

carrier  wavelength 

N(b,B ) 

Gaussian  process  with  the  mean  vector  b  and  covariance  matrix  B 

n 

overall  transition  matrix  of  sight  states 

n, 

transition  matrix  of  the  /'- th  AN’s  sight  state 

( Qh ) 

range  offset  variance  (covariance  matrix) 

% 

carrier  phase  of  AN  transmitter  at  time  of  transmission 

cpit) 

partial  carrier  phase  cycle  measurement 

x1 

chi-square  distribution 

dim(.) 

dimension  function 

Re{.} 

real  part 

V. 

Jacobian  function 
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Chapter  1  Overview  of  contract  objective,  phases  and  deliverables 

1.1  Overview  of  the  developed  solution  of  SGL 

The  overall  objective  of  the  research  activities  comprising  this  contract  is  mainly  that  of  the 
development  of  self  geo-location  (SGL)  in  which  the  mobile  unit  (denoted  as  the  UE  for  user 
equipment)  is  to  be  located  relative  to  the  surrounding  environment  based  on  using  LTE 
wireless  signals  of  opportunity.  To  facilitate  this  we  consider  inputs  as  shown  in  Figure  1-1. 


prior  information  and  maps 


Signals  from  LTE  and  other  wireless  sources 


GNSS  signals  when  available 


Camera  sensor 


Other  sensors 

-  accelerometers,  rate  gyros,  magnetometers 


Figure  1-1  Inputs  to  the  UE  for  SGL 


Recent  technology  articles  have  outlined  the  race  to  have  "Google  Maps"  extended  to  indoor 
locations,  malls,  airports,  arenas  and  so  forth.  Smartphone  manufactures  such  as  Apple,  Nokia, 
Samsung  and  RIM  and  ASIC  manufactures  such  as  Infineon,  Qualcomm  and  Broadcom  are 
scrambling  to  achieve  processing  that  will  provide  sufficient  resolution  such  that  the  detail  of 
these  available  maps  will  be  useful  to  the  smartphone  user.  Such  technology  has  obvious 
application  for  military  based  personal  and  asset  location.  Differences  will  be  that  the  wireless 
signals  of  opportunity  that  can  be  used  for  such  positioning  are  perhaps  uncooperative  with 
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poorly  defined  source  locations  and  no  guarantee  of  being  accurately  synchronized.  The  overall 
objective  of  the  research  is  the  enhancement  of  signal  processing  algorithms  for  estimating  the 
location  of  a  mobile  station  (MS)  or  user  equipment  (UE)  based  on  wireless  signalling.  The 
wireless  signalling  is  assumed  to  be  sourced  from  4G  Long  Term  Integration  (LTE)  access  nodes 
(AN)  which  are  primarily  used  for  mobile  data  and  voice  communications.  The  downlink 
signalling  is  specifically  exploited,  as  the  desire  is  for  the  UE  processing  to  be  capable  of  stand 
alone  "self  geo-location"  (SGL).  This  implies  that  the  network  location  processing  or  assistance 
is  not  involved.  Further  it  is  not  assumed  that  the  UE  be  registered  with  the  LTE  network  but 
rather  that  it  exploits  LTE  downlink  signals  of  opportunity  for  the  location  estimate.  Hence  the 
relevant  signalling  will  come  from  the  LTE  downlink  reference  signals  (RS)  that  are  synchronised 
across  the  LTE  network  of  AN's.  The  signal  structure  of  the  RS's  are  assumed  to  be  known  by 
the  UE  as  well  as  the  location  of  the  AN's. 

Additional  inputs  can  come  from  GNSS  (GPS)  signals  when  they  are  available.  The  camera  input 
from  a  small  web-cam  like  device  integrated  into  the  UE  has  been  shown  to  be  of  very  valuable 
input.  Other  sensors  are  the  accelerometers,  gyros  and  magnetometers  which  provide  some 
information  relevant  to  SGL.  In  this  work  development  effort  has  been  expended  on  the 
ancillary  inputs  provided  by  the  camera  and  accelerometers. 

The  simplest  scenario  for  the  UE  SGL  within  an  LTE  signalling  environment  is  where  the  ANs  are 
in  known  locations  and  synchronized  in  terms  of  the  carrier  frequency  and  code  phase.  This  is 
typical  of  the  CDMA  IS2000  network  that  was  studied  in  a  previous  contract  in  which  the  code 
phase  synchronization  tolerance  is  generally  better  than  50  nsec.  The  code  offset  of  the  pilot 
channel  actually  identified  the  particular  base  station,  hence  the  code  phase  had  to  be 
synchronized  to  facilitate  this.  Also  the  IS2000  transmissions  were  calibrated  in  terms  of  delay 
such  that  they  were  within  a  50  nsec  tolerance  relative  to  GPS  time.  This  facilitated  network 
based  location  services  with  the  initial  impetus  being  compliance  with  the  e911  requirements. 
However,  this  is  an  added  complexity  that  needs  to  be  maintained  that  LTE  dispenses  with. 
Clearly  in  LOS  propagation  environments,  the  SGL  tracking  is  robust  provided  that  there  are  a 
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sufficient  number  of  ANs  available  and  that  the  geometry  does  not  significantly  compromise  the 
GDOP  as  is  illustrated  in  Figure  1-2.  Here  the  UE  is  frequency  and  time  synchronized  with  the 
set  of  AN’s  that  are  all  located  such  that  there  is  a  LOS  connection  with  the  UE.  Furthermore 
the  AN  locations  are  precisely  known  to  the  UE.  It  is  assumed  that  the  AN  to  UE  ranges  can  be 
determined  accurately.  For  the  particular  scenario  depicted  in  Figure  1-2  there  are  only  two 

unknowns,  namely  the  (x,_y)  position  of  the  UE  and  there  are  4  range  measurements  possible 

from  the  synchronized  AN’s  to  the  UE.  Hence  a  single  isolated  position  estimate  can  be 
achieved  with  redundancy.  Not  prior  trajectory  knowledge  or  assumptions  of  the  UE  is 
necessary. 


Figure  1-2  Simplest  case  of  SGL  with  all  synchronized  AN's  in  known  locations  and  LOS 
propagation 

The  implied  assumption  that  the  UE  and  AN  clocks  are  synchronized  is  of  course  invalid  as 
these  clocks  are  independent  with  no  direct  physical  or  implied  connection.  Hence  it  is 

necessary  to  add  the  relative  UE  clock  bias  to  the  variable  list  such  that  we  have  ( x,y,b )  as  the 

variable  list.  This  type  of  SGL  is  equivalent  to  the  GNSS  navigation  solution  or  the  SGL  within  a 
terrestrial  CDMA  IS2000  network.  Note  that  in  the  scenario  depicted  we  still  have  a 
measurement  redundancy  allowing  for  an  isolated  position  estimate  of  the  UE. 
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Tight  timing  synchronization  amongst  the  LTE  AN’s  is  not  necessary  for  optimum  data 
communication  functions  nor  is  it  necessary  for  facilitating  network  based  location,  hence  there 
is  little  impetus  for  the  service  provided  to  accurately  time  synchronize  the  network. 
Consequently,  assuming  that  the  bias  of  the  code  phase  of  the  AN  emissions  is  sufficiently 
accurately  synchronized  for  SGL  is  not  valid.  Hence  it  is  necessary  to  carry  a  state  variable  for 
each  observable  AN  for  the  code  phase  bias.  However,  it  is  a  reasonable  assumption  that  the 
LTE  AN  signals  are  frequency  synchronized  as  such  reference  signals  are  available  from  GNSS 
sources  and  back-haul  links.  As  determined  in  phase  I  simulations,  allowing  for  a  static  bias  of 
the  code  phase  of  each  AN  did  not  significantly  degrade  the  robustness  of  the  SGL  tracking 
even  though  there  were  more  state  variables  to  contend  with.  However,  we  now  have  more 
variables  than  observables  for  a  single  static  location  of  the  UE.  A  further  complexity  is  that  it  is 
not  guaranteed  that  the  LTE  AN  is  in  a  precisely  known  location  which  is  likely  the  case  in  a 
non-cooperative  scenario.  This  issue  in  addition  to  not  knowing  the  code  phase  implies  that 
isolated  UE  location  estimates  are  not  possible. 

This  significantly  complicates  SGL  in  that  we  have  to  generate  additional  measurements  to 
exceed  the  number  of  unknowns.  A  key  is  that  the  AN’s  are  assumed  to  be  in  static  positions 
and  are  frequency  locked.  Consequently  if  two  consecutive  measurements  of  the  AN-UE  links 
are  taken  then  we  can  amass  more  measurements  than  unknowns.  However,  the  consecutive 
sets  of  AN-UE  measurements  are  clearly  linearly  dependent  and  hence  to  not  lead  to  a 
sufficient  set  of  linearly  independent  measurements  necessary  for  resolving  all  of  the  unknowns. 
However,  if  the  AN  is  moved  between  the  sets  of  UE-AN  measurement  sets  then  independent 
measurement  sets  result.  This  is  illustrated  in  Figure  1-3. 
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Figure  1-3  Multiple  sets  of  measurements  generated  by  moving  the  UE 

Assume  that  a  set  of  K  UE  positions  are  used  and  that  for  each  position  a  complete  set  of 
measurements  is  available.  If  each  measurement  set  is  entirely  independent  then  there  is 
nothing  gained  from  combining  the  K  sets  of  measurements.  However,  the  trajectory  of  the  UE 
is  typically  quite  smooth  and  also  the  random  instability  of  the  AN  and  UE  clocks  are  highly 
correlated.  Hence  while  we  have  not  explicitly  resolved  the  issue  that  the  number  of  unknowns 
exceeds  the  number  of  measurements,  then  at  least  we  have  minimized  the  uncertainty  of  the 
variables.  Consider  it  this  way.  If  noise  free  measurements  are  available  that  exceed  the 
number  unknowns  then  deterministic  calculations  of  the  unknown  UE  position  variables  are 
possible.  However,  the  measurements  have  a  random  component  hence  the  estimate  of  the 
UE  position  will  be  random.  With  taking  the  K  strongly  correlated  measurements  into  account 
where  the  UE  is  moved  into  K  positions  along  a  smooth  trajectory,  the  issue  is  to  exploit  all  of 
the  correlations  that  are  there  as  deterministic.  The  end  result  is  that  the  moving  UE  just 
makes  it  more  complicated  to  extract  and  exploit  all  of  the  underlying  correlations. 

Consequently  we  need  a  method  that  could  efficiently  use  any  prior  information  and  disparate 
measurements  that  were  related  to  the  unknown  variables  in  arbitrary  ways.  The  only  means 
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of  robust  SGL  in  such  a  context  is  to  use  Bayesian  filtering  methods.  The  Kalman  filter  (KF), 
extended  Kalman  filter  (EKF)  and  the  particle  filter  (PF)  have  been  thoroughly  investigated  in 
this  context  as  part  of  this  research.  The  complexity  of  the  problem  revealed  that  additional  tools 
were  needed  to  accommodate  the  numerous  and  disparate  constraints  and  observations.  This 
was  really  the  impetus  for  involving  the  Simultaneous  Location  and  Mapping  (SLAM)  algorithms 
that  were  developed  for  the  LTE  SGL.  It  was  discovered  that  in  many  scenarios  of  practical 
interest  that  the  UE  could  use  signals  of  opportunity  from  the  LTE  sources.  That  is,  just 
knowing  that  the  LTE  AN  transmission  comes  from  a  stationary  source  is  sufficient  to  gainfully 
use  it  in  a  SLAM  context.  Even  if  all  of  the  AN’s  are  in  unknown  locations  then,  in  principle,  it  is 
possible  to  simultaneously  locate  the  AN’s  and  the  trajectory  of  the  UE.  Also  a  more  efficient 
fastSLAM  method  was  developed  for  the  SGL.  An  overview  of  these  Bayesian  methods  is 
given  in  chapter  2  with  the  details  presented  in  part  II  of  this  report. 

To  proceed  with  a  Bayesian  approach  we  make  the  justifiable  assumption  that  the  AN’s  are 
assumed  to  be  stationary  and  that  the  AN  signals  are  frequency  locked.  Note  that  this  does  not 
imply  that  they  are  time  synchronized  such  that  the  code  phase  is  arbitrary  but  the  offsets  are 
relatively  static  in  comparison  to  the  AN  emissions.  A  separate  bias  for  each  AN  emission  is 

then  assumed  and  denoted  by  the  variables  which  are  unknown  but  constant  with 

time.  In  Figure  1-3,  the  UE  is  moved  while  taking  range  measurements  with  the  underlying 
assumption  that  a  Bayesian  estimator  can  be  realized  that  incorporates  the  assumption  of  a 
smooth  trajectory  with  an  uncertainty  imposed  at  each  time  step.  The  Bayesian  filter  can  also 
take  into  account  the  prior  knowledge  of  the  UE  position  from  the  last  known  position  or  user 
input  relative  to  a  supplied  map. 

A  significant  reduction  in  the  uncertainty  of  each  of  the  UE  update  steps  is  possible  if  computer 
vision  (CV)  inputs  are  used  from  a  small  webcam  that  can  easily  be  integrated  into  the  handset 
device.  As  described  in  chapter  3,  CV  provides  an  accurate  measurement  of  short  trajectory 
lengths  on  the  order  of  several  tens  of  meters.  Figure  1-4  outlines  the  basic  concept  of  using 
CV  ego-motion  relative  to  a  ground  surface  for  generating  updates  for  the  Bayesian  filter.  Even 
if  these  inputs  drift  over  longer  trajectory  lengths  they  significantly  reduce  the  update  uncertainty 
that  the  Bayesian  filter  has  to  otherwise  contend  with.  Note  that  the  CV  measurements  add 
tighter  constraints  between  the  UE  positions  by  removing  much  of  the  uncertainty  of  the  UE 
update. 
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Figure  1-4  Ego-motion  of  camera  sensor  relative  to  a  2D  manifold  surface 


A  challenge  is  that  the  CV  problem  is  not  trivial  as  the  UE  is  hand  held  and  hence  is  subject  to 
random  motion  in  6  DOF  (3  translation  and  3  orientation  angle  variables).  If  it  can  be  assumed 
that  a  flat  levelled  floor  is  being  observed  then  the  two  tilt  angles  of  the  camera  can  be 
measured  absolutely.  As  such  only  the  azimuth  angle  and  the  height  of  the  camera  needs  to  be 
added  as  nuisance  variables  to  the  state  variable  list.  To  robustly  measure  these  tilt  angles  and 
hence  provide  a  direct  CV  based  trajectory  ego-motion  output  that  can  be  used  in  the  Bayesian 
estimator,  a  Kinect  type  of  projector/camera  combination  is  considered.  A  stereoscopic  camera 
was  also  considered  and  partially  developed  for  this  purpose.  These  means  were  demonstrated 
to  work  extremely  well  however,  the  additional  hardware  complexity  is  a  factor.  As  well  it  does 
not  fit  the  handheld  form  factor.  Hence  monocular  webcam  solutions  are  currently  being 
considered  as  no  additional  hardware  is  required. 

Additionally  3D  CV  perspective  transformations  have  been  developed  in  phase  II  as  providing  a 
potential  solution  to  resolve  the  unknown  nuisance  variables  associated  with  the  UE/camera 
trajectory.  A  promising  possibility  for  absolute  trajectory  start  point  referencing  is  that  of  the  CV 
recognizing  a  wall  mounted  icon,  template  or  other  geometric  structure.  The  perspective 
transformation  is  used  in  this  case  as  is  illustrated  in  Figure  1-5.  The  2D  template  is  shown  on 
the  right  side  with  the  perspective  on  the  left  side.  From  the  vantage  point  of  the  camera  the 
perspective  transformation  can  be  determined  from  which  the  6  DOF  variables  of  the  relative 
camera  position  can  be  determined. 
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camera 


Figure  1-5  Inverse  perspective  transformation  obtained  of  a  known  template  on  the  wall 

CV  processing  with  the  standard  web  camera,  which  is  compatible  with  a  handset  device,  is 
sufficiently  capable  of  estimating  this  perspective  transformation.  Naturally,  if  such  templates  or 
icons  are  readily  observable  throughout  the  UE  trajectory  then  one  can  question  the 
requirement  for  the  LTE  signalling.  However,  icons  are  not  always  available,  in  known  locations 
or  sufficiently  distributed  that  they  are  visible  throughout  the  UE  trajectory.  Certainly  in  airport 
lobbies  or  commercial  buildings  such  icon  systems  will  inevitably  be  implemented.  However,  in 
more  hostile  or  unknown  territory  scenarios,  the  existence  of  usable  icons  is  not  guaranteed  but 
is  probable.  Hence  the  development  has  been  extended  to  markers  of  opportunity  as  observed 
in  the  FOV  of  the  camera.  Also  as  will  become  evident  is  that  the  CV  observables  are  useful  in 
obtaining  accurate  assessment  of  short  trajectory  lengths  but  that  ego-motion  based  on  CV 
alone  drifts.  Hence  the  combination  of  CV  and  LTE  wireless  observables  are  complementary. 

A  major  limitation  with  SGL  based  on  LTE  signals  is  the  multipath  propagation  effects  which 
typically  randomly  disperses  the  signal  such  that  delay  spread  can  be  up  to  several  hundred 
nsec.  Hence  the  deviation  of  the  range  measurements  can  be  several  10’s  of  meters  due  to 
multipath  alone.  This  renders  the  LTE  signals  of  very  little  information  content  in  the  SGL 
context.  Certainly  the  wide  20  MHz  bandwidth  helps  as  the  resolvable  components  of  the 


17/143 


multipath  can  be  eliminated  leaving  only  the  initial  group  of  multipath  components 
corresponding  to  the  leading  edge  and  hence  reducing  the  range  uncertainty  to  about  30 
meters.  Also  there  is  progress  in  LTE  technology  development  to  push  the  bandwidth  up  to  100 
MHz.  However,  with  current  bandwidths,  the  delay  spread  in  NLOS  environments  is  such  that 
the  uncertainty  in  the  pseudo-range  measurement  is  too  large  to  be  of  practical  use  for  indoor 
SGL. 

In  phase  II  consideration  was  given  to  representing  the  multipath  delay  bias  with  a  random 
process  as  the  bias  will  change  between  LOS  and  NLOS  conditions.  The  alternative  is  a 
Boolean  switch  that  indicates  states  of  high  LOS  and  high  NLOS.  The  detector  processing  for 
such  a  switch  is  based  on  the  synthetic  array.  Note  that  the  AN-UE  propagation  path  will  be 
LOS  in  certain  locations.  The  SA  can  be  used  to  accurately  access  the  spatial  coherency  of  the 
incident  LTE  wireless  signal  as  illustrated  in  Figure  1-6.  Here  the  spatial  coherency  of  the  AN 
signal  is  evaluated  as  the  UE  is  moved.  If  the  spatial  coherence  metric,  denoted  byyLOS ,  is  high 
then  a  LOS  condition  is  declared. 
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Figure  1-6  Illustration  of  the  boundary  between  the  NLOS  and  the  LOS  region. 

Based  on  this  assessment  various  options  are  available. 

1 .  Trigger  a  direct  calculation  of  the  estimate  of  the  AN  location  relative  to  the  UE  trajectory  if 
yLOS  is  sufficiently  high. 

2.  If  a  radio  map  exists  that  is  based  on  /£OSthen  modify  the  belief  map  of  the  location  of  the  UE 

in  a  Bayesian  sense  that  takes  the  overall  UE  trajectory  into  account. 

3.  It  is  necessary  to  use  a  state  variable  in  the  SGL  SLAM  that  is  related  to  the  delay  bias 
between  the  UE  and  the  AN.  Different  values  of  the  delay  bias  state  variable  can  be  used 
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depending  on  the  condition  of  multipath.  Hence  yws  can  be  used  to  select  the  appropriate 

model  of  the  delay  bias  for  the  currently  encountered  propagation  condition. 

4.  Multiple  pseudo  range  measurements  are  made  at  each  location  of  the  UE  as  it  moves  along 
the  trajectory.  The  quality  of  these  measurement  in  terms  of  the  SGL  objective  is  generally 
higher  for  LOS  vs  NLOS  conditions.  Hence  yws  can  be  used  to  form  a  weighting  factor  related 
to  how  reliable  a  given  UE-AN  pseudo-range  measurement  is. 

5.  RSS  measurements  are  easier  to  generate  by  the  UE  and  can  also  be  incorporated  as  a  set 
of  measurements  for  the  SGL.  If  yws  is  low  then  a  predominantly  NLOS  condition  exists  in 

which  the  RSS  observables  often  provide  more  information  than  the  pseudo-range 
measurements.  Hence  the  weighting  of  point  4  can  be  extended  to  the  RSS  observables. 

Figure  1-7  illustrates  the  trajectory  of  the  UE  intercepting  the  spherical  wave  emanating  from  the 
AN.  If  the  local  trajectory  can  be  accurately  determined  say  from  a  CV  input  then  if  yLOS  is  such 

that  there  is  a  LOS  condition  then  it  is  possible  to  determine  the  location  of  the  AN  relative  to  the 
trajectory  as  stated  resulting  in  a  high  information  observable. 


Figure  1-7  Assessment  of  the  spatial  coherence  by  continuous  sampling  of  the  AN  complex 
gain 
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An  overall  flow  chart  of  the  SGL  processing  is  given  in  the  Figure  1-8.  Initially  prior  maps  and 
information  are  used  by  an  initialization  phase  which  configures  the  SLAM  algorithm  for  SGL. 
The  initial  displacement  increment  is  taken  along  the  trajectory  and  the  observables  are  sampled. 


Figure  1-8  Overall  flow  of  SLAM  algorithm 

Initialization  phase  -  the  UE  is  enabled  for  the  SGL  operation.  Any  prior  maps  are  read  into 
the  device  along  with  knowledge  of  where  it  was  last  time  it  did  an  SGL  estimate.  Any  user 
input  information  is  of  course  also  applied.  Also  radio  maps  are  loaded  in  corresponding  to  the 
approximate  location.  These  will  include  the  list  of  possible  AN’s  in  the  vicinity  that  are  known. 
Additionally,  the  statistical  trajectory  information  will  be  input  and  processed.  This  may  include 
knowledge  of  whether  the  LIE  is  mounted  to  push  cart,  vehicle,  foot  soldier  or  other  asset.  A 
search  of  GPS  and  perhaps  other  GNSS  sources  is  done  with  a  listing  of  the  available  sources  for 
clock  synchronization  and  possibly  AN  sources.  Next  the  UE  clock  offset  is  estimated  from  the 
set  of  observable  signals.  It  is  assumed  here  that  the  UE  is  relatively  stationary  and  that  the 
signals  are  categorized  for  their  quality  in  terms  of  use  as  frequency  reference  signals.  A  model 
for  the  clock  oscillator  instability  is  then  determined  as  will  be  used  in  the  SLAM  method.  Then 
the  UE  camera  is  used  to  scan  for  a  LandMark  (LM)  that  can  correspond  with  the  map.  The  user 
can  optionally  assist  in  this  stage  by  scanning  the  camera.  Of  course  if  there  are  no  relevant 


20/143 


maps  available  then  this  step  is  skipped.  Based  on  the  observed  CV  LM’s  and  the  set  of  GNSS 
signals  and  AN  signals  that  are  available  the  SLAM  algorithm  is  configured  along  with  an  initial 
belief  map  of  the  UE  (x,y)  (and  possibly  z)  position  that  is  used  for  the  SLAM. 


Figure  1-9  Initialization  phase 

Sample  observables 

The  various  observables  are  sampled  as  configured  from  the  set  in  the  initialization  step.  These 
involve  the  wireless  signals  from  the  LTE,  WiFi,  RFID  or  other  sources  of  opportunity,  any 
available  GNSS  signals,  camera  inputs  as  well  as  inputs  from  MEMS  inertial  sensors  and 
magnetometers. 

SLAM  update 

The  fastSLAM  method  is  used  for  the  Bayesian  filter  for  updating  the  trajectory  and  the  feature 
point  information.  The  most  efficient  combination  is  the  PF  used  for  the  UE  displacement  and 
the  EKF  for  tracking  nuisance  variables  and  converging  map  FP’s. 
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Figure  1-10  SLAM  update 

Finally,  the  possibility  of  a  map  is  being  considered  for  the  SGL  activity.  Such  a  map,  if 
available,  would  provide  the  UE  with  an  outline  of  the  building  from  which  a  prior  belief  map  of 
location  can  be  assessed.  It  would  show  the  position  of  permanent  obstacles  and  walls  such  that 
probabilities  of  estimated  trajectories  can  be  assessed.  This  is  a  very  useful  input  for  the  various 
Bayesian  filters  that  are  being  developed.  The  map  could  also  provide  location  of  templates  that 
can  be  observed  by  CV  and  perhaps  the  locations  of  AN’s.  In  hostile  scenarios  where  a  partial 
map  is  available  at  best,  an  outcome  of  the  SLAM  could  be  to  slowly  anneal  such  a  map  and  add 
features  to  it.  This  is  a  powerful  tool  especially  if  there  are  several  UE’s  in  the  area 
cooperatively  building  up  such  a  map. 


1.2  Outline  of  the  report 

The  organization  of  the  remainder  of  this  final  report  is  as  follows: 


22/143 


Chapter  2.  Overview  of  the  SLAM  algorithm  for  SGL  -  This  chapter  describes  the  SLAM 
algorithm  as  applicable  to  the  SGL  problem.  The  general  set  of  observables  from  the  AN 
wireless  downlink  transmissions  and  the  CV  observables  will  be  described  in  the  context  of 
being  inputs  for  the  SLAM  algorithm.  The  PF  from  the  phase  I  will  be  reviewed  followed  by  a 
development  of  the  SLAM  based  on  the  Kalman  filter  (KF),  extended  Kalman  filter  (EKF)  and 
then  the  fastSLAM  which  involves  a  combination  of  the  EKF  and  PF.  fastSLAM  will  be  what  is 
used  specifically  for  the  SGL  algorithm  development.  The  details  of  the  SLAM  algorithm 
research  are  contained  in  Part  B.  Also  in  this  section  is  a  detailed  description  of  the  downlink 
LTE  signal.  As  described,  the  resource  blocks  in  the  reference  channel  contain  the  coherent 
synchronization  blocks  from  which  the  signal  phasing  at  the  mobile  can  be  estimated. 

Chapter  3.  CV  Ego-motion  -  There  have  been  many  new  developments  in  the  CV  ego-motion 
that  are  of  significance  to  the  objectives  of  this  contract.  These  have  been  developed  to  the 
point  that  integration  with  the  SLAM  algorithm  is  practical.  The  CV  forms  an  independent  set  of 
observables  that  is  useful  for  short  term  trajectory  estimation.  Not  only  is  the  CV  ego-motion 
used  as  SLAM  observables  but  also  as  a  means  of  testing  the  algorithms  for  the  LOS  AN  wireless 
beamforming.  Given  a  good  set  of  feature  points,  the  ego-motion  is  sufficiently  accurate  to 
precisely  track  the  translation  and  orientation  of  the  wireless  receiver. 

Chapter  4:  Development  of  Experimental  Hardware  -  Several  iterations  of  the  hardware  have 
been  undertaken  with  the  current  unit,  capable  of  some  simultaneous  recording  and  fusion  of 
wireless  and  CV  signals  in  the  SLAM  context.  The  receiver  processing  consists  of  two  parts.  The 
first  is  the  signal  is  down-converted,  sampled  and  stored.  Then  post-processing  simulates  the 
SLAM  tracking  algorithms  along  with  the  synthetic  beam  forming. 

Chapter  5:  Discussion  -  Overall  observations  and  outcomes  of  the  theoretical  and  experimental 
outcomes  and  suggested  next  steps. 
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Chapter  2  SLAM  SGL  Algorithm 


2.1  Overview 

In  this  chapter  the  outline  of  the  SLAM  method  developed  for  the  SGL  system  will  be  described. 
The  bulk  of  the  details  are  in  part  B  of  the  report  which  describes  the  OWLS  (Opportunistic 
wireless  location  system).  As  described  in  Chapter  1,  the  SGL  problem  involves  the  estimation 
of  the  UE  trajectory  with  LTE  wireless  signals  emanating  from  AN's  from  which  range 
measurements  are  extractable.  AN's  occur  in  two  types:  Anchor  nodes  where  the  location  of 
the  AN  is  known  which  is  denoted  as  ANan  and  feature  nodes  of  initially  unknown  position 
denoted  as  ANfn.  In  either  case  there  are  additional  unknowns  such  as  the  delay  bias 
attributed  to  each  UE-AN  link  as  well  as  the  orientation  of  the  UE  which  affects  the  antenna 
response.  SLAM  processing  is  ideally  suited  for  this  type  of  SGL  problem  and  will  be  developed 
in  this  chapter. 

The  basic  assumptions  are  as  follows: 

1.  LTE  wireless  signalling  will  be  used  for  the  observables 

2.  Only  downlink  LTE  signals  will  be  assumed.  The  UE  is  not  registered  with  the  network  and 
hence  uses  the  LTE  downlink  signals  as  signals  of  opportunity. 

3.  The  LTE  downlink  signals  from  the  various  AN's  may  or  may  not  be  accurately  time 
synchronized.  Both  modes  will  be  considered. 

4.  The  UE  has  to  have  all  the  information  regarding  the  SGL  estimate  such  that  a  handset  based 
location  solution  is  possible  for  an  uncooperative  network. 

5.  Multipath  will  be  variable  as  indoor  and  dense  development  areas  are  of  primary  interest. 

6.  The  UE  is  expected  to  be  undergoing  some  arbitrary  motion  that  varies  in  terms  of 
smoothness. 

7.  UE  is  not  frequency  synchronized  to  the  AN  in  any  way.  That  is  it  must  be  assumed  that  the 
UE  has  a  free  running  clock. 

8.  CV  inputs  are  available  for  updates  of  the  UE  pose. 
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9.  SGL  algorithm  must  accommodate  the  possibility  that  the  AN  locations  are  not  accurately 
known. 

The  overall  objective  is  to  generate  a  belief  map  of  the  location  of  the  UE  as  the  trajectory 
unfolds.  The  belief  map  will  provide  a  posterior  estimate  of  the  location  of  the  UE  that  is 
iterated  every  time  interval  which  is  based  on  all  available  data  up  to  the  current  time.  For 
such  estimation  problems  based  on  the  aforementioned  assumptions,  it  is  generally  not  valid  to 
assume  Gaussian  statistics.  Consequently  Kalman  filter  (KF)  type  parametric  Bayesian  filters 
(BF)  are  typically  sub-optimal.  The  particle  filter  (PF)  has  been  developed  for  the  SGL  problem 
as  described  in  the  first  interim  report.  However,  a  limitation  is  that  the  number  of  state 
variables  is  large  such  that  a  mix  of  EKF  and  PF  is  necessary  which  is  the  fastSLAM  algorithm 
that  will  be  developed  for  the  SGL  problem. 

Indoor  environments  are  subject  to  severe  multipath  which  significantly  distorts  wireless 
signals  in  the  LTE  frequency  bands.  Even  though  LTE  has  a  potentially  large  20  MHz  bandwidth, 
this  is  typically  not  sufficient  to  provide  adequate  accuracy  of  time  of  flight  measurements  for 
indoor  SGL.  Multipath  distortions  it  itself  renders  traditional  isolated  location  estimate 
methods  of  SGL  virtually  useless  for  indoor  applications  where  an  accuracy  of  1  meter  is 
required.  Combatting  the  effects  of  multipath  is  the  most  significant  challenge  in  providing 
adequately  accurate  LTE  based  SGL.  This  is  exasperated  by  the  fact  that  the  AN's  are  not 
necessarily  synchronized  nor  in  accurately  known  positions.  This  implies  that  there  is  a  random 
delay  offset  that  must  be  considered  with  regards  to  each  AN  that  cannot  be  considered  to  be  a 
Gaussian  random  variable.  The  only  way  to  mitigate  these  limitations  is  if  the  relatively  smooth 
trajectory  of  the  UE  is  exploited  by  the  BF.  With  a  smooth  trajectory,  multiple  measurements 
between  the  UE  and  the  same  AN  can  be  combined.  The  other  key  point  is  that  the  AN  is 
assumed  to  be  stationary.  Additionally  it  is  optionally  assumed  that  CV  inputs  exist  such  that 
the  relative  short-term  update  motion  in  terms  of  translation  and  orientation  can  be 
determined  and  passed  to  the  BF  as  statistical  UE  update  information.  The  CV  will  be  assumed 
to  consist  of  a  small  webcam  which  is  compatible  with  the  monocular  3D  vision  systems.  An 
overall  imposed  constraint  is  that  the  UE  needs  to  be  self-contained  and  that  the  apparatus  and 
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processing  required  for  the  SGL  function  must  be  commensurate  with  the  handheld  form  factor 
hence  the  monocular  limitation. 


Initially  pure  LOS  links  will  be  assumed  between  the  UE  and  the  various  AN's.  However,  these 
will  be  supplanted  with  NLOS  multipath  links  which  will  be  accounted  for  by  assuming  a  time  of 
flight  (TOF)  bias  that  varies  with  time.  Hence  it  is  not  really  a  bias  but  more  a  slowly  varying 
random  variable  representing  the  bias.  This  bias  can  switch  quickly  when  moving  from  LOS  to 
NLOS  conditions.  As  described  previously,  Boolean  switch  models  and  Markov  models  have 
been  analysed  for  the  fluctuating  delay  bias  issue.  We  have  also  analysed  the  possibility  of 
generating  a  metric  based  on  the  equivalent  of  the  Ricean  K  factor  which  provides  an 
approximate  measure  of  the  level  of  multipath  in  the  LTE  wireless  signal  sources. 

Boolean  multipath  model  switches  are  compatible  with  the  Bayesian  filter  however  with  the 
expense  of  carrying  an  additional  state  variable.  An  example  of  a  scenario  where  the 
propagation  changes  from  LOS  to  NLOS  is  shown  in  Figure  2.1.  The  concept  would  be  that  the 
processing  would  have  to  recognize  that  there  has  been  a  change  in  the  propagation  perhaps 
from  the  coherency  of  the  carrier  phase  as  the  UE  propagates.  This  is  where  the  CV  updates 
will  have  a  significant  impact  as  they  will  provide  a  short  term  trajectory  that  is  sufficiently 
accurate  to  facilitate  beamforming  on  the  AN  signal  to  determine  it's  spatial  coherency.  Hence 
the  quality  of  the  LOS  component  can  be  determined  and  different  models  can  be  switched  into 
the  SGL  SLAM  algorithm. 
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Figure  2.1  Illustration  of  the  boundary  between  the  NLOS  and  the  LOS  region.  In  theory,  as  the 
UE  moved  across  this  boundary,  the  bias  model  would  switch  from  LOS  to  NLOS  with 
appropriately  selected  delays 

The  LOS/NLOS  state  variable  can  take  several  forms.  The  first  is  a  Boolean  indicator  whether 
the  signal  carrier  phase  is  sufficiently  coherent  such  that  a  LOS  threshold  is  surpassed.  Hence 
the  Bayesisan  filter  can  incorporate  both  the  LOS  and  NLOS  propagation  models  for  the  signal. 
For  the  LOS  case  a  fixed  excess  propagation  delay  can  be  assumed.  In  the  pure  NLOS  case,  any 
time  delay  observables  are  of  little  value  to  the  overall  SLAM  and  the  RSS  observables  are  used 
instead.  While  individual  RSS  measurements  are  low  in  information  value,  the  sequence  of  RSS 
measurements  can  give  approximate  location  information  where  TOA  based  observables  are 
useless. 

Another  implementation  is  that  the  excess  propagation  delay  is  a  state  variable.  The  receiver 
demodulates  the  resource  blocks  of  the  LTE  signal  which  are  narrow  bandwidth  coded  signal 
bursts  that  are  then  Fourier  transformed  to  obtain  the  channel  time  domain  impulse  response. 
The  TOA  estimate  of  a  given  LTE  signal  is  then  taken  as  the  time  of  the  leading  edge  of  the 
transformed  impulse  response  offset  by  the  parameter  "b"  which  is  determined  by  the  PF.  The 
higher  the  aggregate  bandwidth  of  the  LTE  downlink  signal  the  sharper  the  leading  edge 
transition  of  the  estimated  channel  impulse  response  will  be  which  will  result  in  less  variation  of 
"b".  In  turn,  this  implies  that  the  spread  of  the  posterior  pdf  of  the  state  variable  belief  is  more 
confined  resulting  in  an  overall  lower  position  variance. 

In  summary,  the  key  to  indoor  positioning  based  primarily  on  wireless  signalling  is  that  the 
there  is  a  sequence  of  correlated  measurements.  That  is  the  trajectory  of  the  UE  is  relatively 
smooth.  This  correlates  the  observables  such  that  even  if  the  trajectory  is  random,  a  steady 
state  covariance  of  the  state  variable  vector  is  achieved. 

The  Bayesian  Fisher  information  matrix  (BFIM)  is  derived  as  a  central  key  to  the  SLAM-based 
opportunistic  wireless  localization  system  (OWLS)  that  unifies  all  information  from  observables, 
the  MN  trajectory,  and  a  priori  knowledge  in  one  single  matrix.  The  BFIM  illustrates  that 
observability  of  a  SLAM-based  OWLS  is  achieved  under  the  assumptions  of  stationary  ANs  and 
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the  MN  smooth  trajectory.  The  development  of  OWLS  is  convoluted  and  detailed  in  places. 
Consequently,  in  an  attempt  at  a  coherent  overview,  the  details  of  OWLS  are  deferred  to  part  B 
of  this  report. 

What  makes  indoor  positioning  difficult  in  the  LTE  SGL  scenario  is  that  there  are  so  many 
unknowns.  If  the  UE  trajectory  was  the  only  unknown  with  the  AN's  synchronized  and  in  known 
locations,  and  the  propagation  LOS  then  the  SGL  would  be  that  of  the  conventional  GPS 
navigation  problem  which  is  very  well  understood  and  generally  robust.  Unfortunately,  in  the 
LTE  SGL  context,  the  location  of  the  AN's  is  generally  unknown  as  is  the  code  phase  of  the 
signal.  Furthermore  the  multipath  propagation  results  in  a  TOA  bias  that  is  unknown  and  is  not 
stationary.  Further  the  clock  of  the  UE  is  not  synchronized  to  any  reference  point  in  the 
network.  Hence  SGL  is  a  problem  where  the  UE  trajectory  is  a  small  component  of  the  overall 
list  of  unknowns  with  the  bulk  of  the  unknowns  being  the  vast  list  of  nuisance  parameters 
involving  the  AN's  and  the  multipath  propagation.  For  this  reason  the  development  of  SLAM 
(Simultaneous  Localization  and  Mapping)  has  been  pursued.  SLAM  is  a  promising  algorithm  for 
this  situation  as  the  mapping  component  can  be  more  abstract  than  merely  determining  the 
position  of  the  AN's  as  part  of  the  map  discovery.  SLAM  will  simultaneously  determine  the 
unknown  parameters  of  the  AN's  and  the  trajectory  of  the  UE  itself.  Hence  the  "map" 
component  consists  of  the  long  list  of  nuisance  parameters.  The  key  result  of  the  early  work  in 
SLAM  was  that  there  is  a  high  degree  of  correlation  between  estimates  of  the  location  of 
different  landmarks  in  a  map  and  these  correlations  would  increase  with  future  observations 
[1].  Later  in  decade  of  1990,  the  key  research  on  convergence  were  developed  by  Csorba  [8] 
and  then  on  KF  based  SLAM  methods  and  the  probabilistic  localisation  and  mapping  methods 
by  Thrun  [9].  The  conceptual  breakthrough  in  SLAM  research  formulated  around  this  time  was 
that  it  was  generally  convergent.  It  was  recognized  that  the  correlations  between  landmarks 
were  the  critical  component  that  enabled  convergence. 
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To  describe  SLAM  in  the  SGL  context  begin  with  a  consideration  of  the  isolated  UE  for  which  we 
wish  to  estimate  the  position.  Clearly  this  is  impossible  without  any  observations  so  assume  an 
AN  that  generates  a  recognizable  pilot  signal  as  observed  by  the  UE  which  perhaps  does  a  TOA 
measurement  of  the  signal  relative  to  it's  own  clock.  Unfortunately  this  provides  no 
information  as  the  AN  location  is  not  known,  the  downlink  signal  is  not  synchronized  and  the 
phase  of  the  UE  clock  is  random.  Now  suppose  that  the  UE  moves  and  at  an  other  time  instant 
observes  the  AN  signal.  As  the  two  UE  positions  are  assumed  to  be  independent,  nothing  is 
gained  from  the  second  observation.  If  the  UE  was  stationary  such  that  the  measurements 
were  totally  correlated,  we  would  just  have  two  measurements  of  the  same  thing.  In  fact  the 
measurements  would  still  be  uncorrelated  as  the  clocks  are  not  synchronized  and  that  the  time 
interval  between  measurements  is  assumed  to  be  entirely  random. 

But  now  suppose  that  the  time  between  measurements  is  a  fixed  time  set  by  the  UE  and  that 
the  frequency  offset  and  drift  between  the  clocks  has  been  noted  by  the  UE.  Then  based  on  the 
sequence  of  two  measurements  it  is  possible  to  say  that  the  UE  radial  distance  from  the  AN  has 
changed  by  an  estimated  amount  subject  to  some  uncertainty  traceable  to  the  relative 
instability  of  the  clocks  and  perhaps  TOF  variations  due  to  multipath.  However,  this  is 
information  regarding  the  UE  trajectory  albeit  modest.  Now  suppose  that  the  UE  is  moved  by 
an  operator  who  is  constrained  to  some  form  of  stochastic  model.  That  is  we  know  the 
operator  is  walking  and  not  in  a  high  speed  vehicle  for  example.  Now  the  radial  position 
information  can  be  combined  with  the  user  motion  dynamics  model  to  present  a  refined 
estimate  regarding  the  trajectory.  Now  consider  this  as  extended  over  a  number  of 
observations  as  the  UE  is  moved.  We  can  infer  plausible  trajectories  from  the  sequence  of 
observations  and  the  prior  assumptions  of  the  dynamics  model.  The  end  result  after  the 
sequence  of  observations  is  a  belief  map  or  PDF  of  the  possible  UE  trajectory  relative  to  the 
unknown  position  of  the  AN.  Even  though  the  absolute  position  is  not  known,  there  is 
quantifiable  information  in  the  resulting  conditional  belief  map. 


29/143 


Next  consider  another  AN  from  which  the  UE  can  make  observations.  The  only  assumption  is 
that  both  AN's  are  stationary.  The  UE  is  moved  while  a  sequence  of  observations  is  made. 
Without  observations,  as  time  progresses,  the  PDF  of  the  UE  location  will  diffuse  due  to 
uncertainties  in  the  UE  motion  itself.  The  observations  themselves  are  subject  to  clock 
instability  and  variations  in  the  multipath  delays.  However,  the  net  information  of  the 
observations  will  tend  to  offset  the  probability  diffusion  that  occurs  due  to  the  UE  dynamics 
itself.  It  is  the  correlation  between  the  sequence  of  measurements  that  provides  the 
information  rather  than  the  isolated  measurements  themselves.  SLAM  will  extract  all  of  the 
information  contained  in  these  correlations  and  establish  a  joint  belief  map  of  all  of  the 
unknown  variables.  These  include  not  only  the  UE  trajectory  but  also  the  nuisance  or  mapping 
variables. 


Figure  2.2  Graphical  schemes  of  a)  online  SLAM  in  comparison  with  b)  full  SLAM 


As  the  UE  is  moved  and  the  location  of  the  M  fp's  and  the  UE  trajectory  is  estimated,  the  errors 
in  the  state  variables  are  correlated  [2].  This  correlation  is  exploited  in  SLAM  and  can  be 
conceptualized  as  an  emerging  mesh  of  constraints  as  illustrated  in  Figure  2.2.  Here  z  denotes 
the  measurement,  m  the  variables  associated  with  the  AN's  and  other  nuisance  parameters,  u 
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denotes  the  UE  position  or  the  trajectory  and  q  the  dynamics  model  of  the  UE  location  update. 
Typically  an  on-line  SLAM  algorithm  is  required  that  continuously  processes  the  information 
and  provides  the  output  of  the  most  current  estimate  of  the  UE  trajectory  and  mapping 
variables.  This  is  illustrated  in  Figure  2.2a.  In  Figure  2.2b  the  off-line  SLAM  is  illustrated  in 
which  the  overall  UE  trajectory  estimate  and  mapping  variable  estimate  is  done  off  line.  Hence 
the  overall  trajectory  is  not  available  until  all  of  the  observations  have  been  made.  While  the 
off-line  gives  better  performance  as  it  uses  all  of  the  measurements  made  it  has  an  obvious 
latency  issue  such  that  it  is  not  suitable  for  real  time  applications.  One  illustration  of  the  SLAM 
constraints  is  that  of  representing  the  individual  constraints  as  springs  between  the  fp's  and  UE 
pose  as  shown  in  Figure  2-2.  An  observation  at  each  time  is  analogous  with  a  displacement  in 
spring  network  and  it  correlation  effects  on  neighbors  depends  on  their  distances  to  other 
feature  points,  i.e.  ,  the  nearer  they  are,  the  greater  is  effect.  As  the  UE  moves  and  received 
more  observations  of  feature  points,  the  correlation  between  feature  points'  estimates 
increase  and  the  spring  become  stiffer.  Noise  and  un-modelled  uncertainty  can  be  imagined  as 
disturbances  to  this  spring  model  that  distorts  the  estimates  of  the  state  variables.  Hence  the 
stiffer  the  springs,  the  more  robust  the  solution  and  the  less  the  influence  of  these  noise 
disturbances  will  cause  errors.  Note  that  the  map  features  could  be  in  known  locations.  As 
such  they  are  anchors  from  which  absolute  location  reference  can  be  inferred.  Also  map 
features  could  have  various  degrees  of  freedom.  For  instance,  a  floor  marker  may  not  be 
unique  but  may  be  oriented  in  a  specific  direction.  A  difficultly  that  is  encountered  with 
features  in  SLAM  is  that  of  establishing  correspondence.  For  active  AN's,  this  is  not  an  issue  as 
the  signal  is  generally  coded  such  that  it  can  be  uniquely  identified  relative  to  the  other  AN 
signals.  However,  if  the  AN  coding  is  not  known  beforehand  then  it  is  difficult  to  separate  it 
from  the  other  AN  emissions  compounding  the  correspondence  problem.  SLAM  maps  for  SGL 
can  extend  beyond  the  AN's  as  CV  is  used.  Hence  while  CV  is  powerful  in  terms  of  establishing 
and  using  map  fp's  there  is  the  issue  of  determining  correspondence  as  will  be  developed  in  the 
following  chapter. 

2.2  SLAM  Information  Matrix 
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SLAM  in  the  SGL  context  is  a  unifying  method  that  can  accommodate  multiple  disparate 
observations  over  multiple  time  intervals.  It  can  also  accommodate  what  is  known  and 
unknown  about  the  various  features  provided  that  the  uncertainty  can  be  quantified  by 
reducing  it  to  a  set  of  covariances.  SLAM  thereby  forms  a  unified  framework  for  fusing  the 
information  from  all  of  the  disparate  bits  of  information,  constraints,  measurements  and  so 
forth  into  a  single  information  matrix.  It  is  therefore  ideal  for  the  LTE/CV/SGL  problem 
formulation. 

The  concept  of  the  information  matrix,  denoted  in  general  by  J ,  is  key  to  the  SLAM  algorithm. 
It  can  loosely  be  regarded  as  a  quantitative  summary  of  all  of  the  information  regarding  the 
unknowns  in  terms  of  the  map  variables  and  the  UE  trajectory  in  one  matrix.  The  relation  to 
quantifiable  information  is  only  valid  for  the  case  where  the  problem  is  linear  and  jointly 
Gaussian  however  for  problems  where  the  nonlinearity  is  mild  and  that  the  measurement  and 
update  noise  is  unimodal  and  approximately  Gaussian  then  the  information  matrix  has  deep 
significance.  To  develop  the  information  matrix  assume  that  the  MLE  of  a  set  of  variables  of  q 

can  be  approximated  as  being  jointly  Gaussian  represented  by  Here  q  represents 

the  actual  value  of  the  state  variable  vector  expressed  as 


<h 


q  = 


q, 


(i) 


_m _ 

For  typical  problems  that  are  not  strictly  Gaussian,  the  approximation  of  the  PDF  of  the  MLE  as 
is  still  generally  valid.  Regardless,  J  is  the  information  matrix  and  is  the  inverse  of 

the  covariance  matrix  of  the  MLE  estimators  for# .  The  concept  of  the  information  matrix  is 
best  exemplified  through  the  following  set  of  simple  examples. 


Consider  a  case  where  the  UE  undergoes  random  motion  in  the  x  direction  according  to  a  first 
order  markov  (random  walk)  model.  Consider  the  initial  step  where  x0  □  0,<Tq  )  which  is  of 

course  the  prior  PDF.  Now  consider  the  step  between  x0  and  x,  as  shown  in  Figure  2-4. 
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Figure  2-4  One  step  problem 


The  update  from  t=0  to  t=1  is  characterized  by 

Xj  =  x0  +  u 
uUN[0,o-2u) 

The  joint  posterior  PDF  of  {x0,Xj}  is  given  as  a  normal  of  zero  mean  and  covariance  of 


Q  = 


2  2  2 

cr0  cr0  +  (7u 


(3) 


Note  that  the  covariance  of  xx  given  x0  is  denoted  as  crL  is  determined  as 


'l|0 


=  (^o  +°l)- 


JL  =  (J2 

_2  u 


which  makes  sense.  Likewise  the  covariance  of  x„ 


x,  is 


(4) 


^-4  ^22 


2  2  2  2 

°o+<r«  o-0+(7u 


(5) 


which  for  the  case  where  c7q  □  a2  is  about  cr2 .  If  there  is  no  measurement  then  the  variance 

of  the  estimate  of  Xj  is  a2  +a2  which  is  also  consistent  with  the  covariance  matrix.  The 
inverse  of  Q  is  given  as 


J  =  Q-1 


2  2 

°0  +  Gu 

^_2  7. 

^0  au 

_ -2 


Note  that  J  can  be  expressed  as 


(6) 
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J  = 


-2  -2  -2 

a  +an  -cr. 


—2 


-2 


-2 


0 


0  0 


-2 


_-2 


-cr 


-2  —2 

-cr  cr,. 


(7) 


which  is  the  combined  information  from  the  initial  condition  and  from  the  update.  Note  that 


[JL  ^olxi 

Consider  another  problem  where  we  have  a  measurement  associated  with  an  anchor  point  of 
ml  as  illustrated  in  Figure  5.  Assume  that  the  variance  of  the  measurement  of  the  anchor 

point  is  crz . 


Figure  2-5  Single  measurement  of  an  initial  trajectory  point 


The  covariance  of  x0  given  the  measurement  is 


_2  _  2  u0 
°0|z  <T> 


cr_crn 


2  2  2  2 

crz+cro  cr_  +  (j0 


(9) 


Now  suppose  that  we  consider  the  measurement  as  an  amendment  to  the  prior  one  step 
problem  as  shown  in  Figure  2-6.  The  two  steps  with  the  measurement  of  the  anchor  feature. 
From  the  graph  we  can  infer 


cr  cr, 


o 


2  2 

+^0 


+  cr; 


(10) 
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Figure  2-6  Single  step  with  a  single  measurement 
Alternately,  consider  the  state  vector  of  q  =  [x0  xjr : 
1 .  first  start  with  the  covariance  matrix  of  the  positions 


0  = 


Gl 

2  2  2 

crn  cr„  +cr 


'o 


'0  w0 

2.  The  information  from  the  measurement  is  given  as 


J  = 


0  0 


as  it  refers  to  the  x0  position  state. 
3.  Compute  the  total  information 


J*=QT1  + 


0 


0  °. 

4.  Determine  the  inverse 


-i 


Qtot  J  tot 

For  the  current  example  we  can  determine  that 


e„  = 


—2  _2 

^0  Gz 


2  2 
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Note  that 


\Qtot\22 


22  22  22 
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(J  z _ z  u _ (J  u  

^2  .  ^2 
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(11) 


(12) 


(13) 


(14) 


(15) 


(16) 
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Note  that  incrementing  the  information  matrix  J  with  a  new  incremental  entity  of  information 
such  as  a  measurement  is  trivial  as  illustrated  by  this  example.  The  only  processing  problem  is 
that  J  has  to  be  eventually  inverted  to  extract  the  actual  variances  of  the  state  variables.  The 
other  issue  is  of  course  that  of  the  underlying  assumption  of  the  problem  being  jointly  Gaussian. 

Note  that  if  the  measurement  provides  no  information  then  cr2  — >  oo  and  cr^0  z  =  cr^  +  cr2  as  in 

the  previous  example. 

Next  suppose  that  independent  measurements  of  the  anchor  were  done  from  both  x0  and  xx 
which  is  illustrated  in  Figure  2-9. 


m 


Figure  2-9  Two  independent  measurements  of  the  position  states  x()  and  Xj 
The  information  matrix  for  the  two  independent  measurements  is  given  as 


Note  that  we  are  implying  here  that  the  noise  or  uncertainty  in  the  measurements  is 
independent  and  hence  the  contribution  to  Jz  is  a  diagonal  matrix.  The  state  covariance  is 
then  determined  as 


(18) 


Next  consider  the  problem  depicted  in  Figure  4-9  but  now  assume  that  the  anchor  point  is  from 


a  feature  point  of  initially  unknown  location.  Now  we  have  a  state  vector  of  q  =  [x0  ^  mf . 

The  covariance  matrix  of  the  updates  is  given  as 
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(19) 


Q  = 


2  2 

<T>  +  *„ 


0 
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The  information  added  by  the  two  independent  measurements  between  m  and  x0  and  between 
m  and  x1  are  captured  in  the  pair  of  information  matrices  as 


Jo  = 


0 


0 

0 


0 


CTz 

0 

1 


(20) 


and 


Ji  = 
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Note  that  while  the  two  measurements  are  independent  they  individually  couple  the  state 
variables  of  m  and  x0  and  m  and  x, .  Hence  the  overall  covariance  matrix  is  given  as 
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An  interesting  observation  in  this  example  is  that  it  is  the  multiple  connectivity  of  a  fp  that  makes 
it  carry  useful  information.  A  fp  measurement  that  is  only  connected  to  one  state  variable  does 
not  convey  information.  We  can  see  this  in  the  present  example  from  the  observation  that 
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Also  it  should  be  pointed  out  that  the  component  of  the  information  matrix  for  the  measurement 
has  the  form  of 


[JArE 

where  p(z  \  q)  is  the  joint  pdf  of  the  measurements.  Note  that  for  the  case  where  the 
measurements  are  jointly  Gaussian  with  a  covariance  that  is  independent  of  q  then 


dqtdq 


-In  p{z\q) 
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Reconsider  the  one  step  problem  shown  in  Figure  4-6  with  the  assumption  that  initially  nothing 
is  known  of  the  position  of  x0  and  .  Now  information  is  given  as  x0  □  N (0,  <j\  )  and 

Xj-XqDA^O,^2).  We  then  have 
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In  the  special  case  where  the  following  conditions  are  met: 

1 .  initial  distribution  of  the  state  variables  is  jointly  normal 

2.  update  of  the  state  vector  is  a  linear  function  of  the  state  at  the  previous  iteration 

3.  measurements  relate  linearly  to  the  state  variables 

then  the  state  variables  conditioned  on  the  measurements  remain  jointly  normal.  In  this  case 
the  MMSE  estimate  of  the  state  variables  given  as  the  conditional  mean  as  (Kay  pg.382) 


qmmse=EU]  +  QM^Z  <26) 

is  the  efficient  estimator.  The  covariance  of  this  estimator  is  denoted  as  Qql_  and  is  given  as 
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(27) 


^q\z  ^qq  ^qz^zz  zZqz 

where  Qq^z  is  the  inverse  of  the  Fisher  Information  matrix.  Hence  the  information  matrix  provides 

a  direct  quantification  of  the  variance  of  the  state  variables  without  having  to  calculate  the 
estimator  first.  Hence  the  information  matrix  provides  an  important  means  of  determining  the 
classes  of  SGL  SLAM  problems  that  are  solvable. 

The  information  matrix  is  the  kernel  of  any  Bayesian  filter  approach  and  especially  SLAM  which 
incorporates  many  disparate  observations.  Recognizing  the  information  associated  with  each 
observation  such  that  it  contributes  constructively  to  the  estimation  of  the  eventual  overall  set  of 
state  variables  is  the  key  to  a  SLAM  implementation.  Details  of  SLAM  based  on  various 
Bayesian  filters  will  be  given  in  part  B. 

2.3  Downlink  LTE  signalling 

In  this  section  we  intend  to  expound  on  the  details  of  the  downlink  of  an  LTE  system.  That  is 
the  link  between  eNodeB  to  UE.  For  convenience  we  also  denote  the  eNodeB  as  the  access 
node  or  AN.  LTE  also  known  as  evolved  universal  terrestrial  radio  access  (E-UTRA),  Super  3G 
(S3G)  has  differentiated  itself  from  WCDMA  by  using  Orthogonal  Frequency  Division 
Multiplexing  (OFDM)  as  a  way  to  transfer  data  for  its  downlink  scheme  and  Single  Carrier  FDMA 
(Frequency  Division  Multiple  Access)  for  its  uplink  scheme. 

LTE  is  designed  to  transfer  packets  and  should  satisfy  the  following  specifications. 


•  100  Mbit/s  for  downlink  and  50Mbit/s  for  uplink 

•  Optimized  operation  for  mobility  of  0  to  15  Km/h  and  high  performance  between 
15Km/h  to  120  Km/h  and  connected  at  speeds  above  120  Km/h  till  350  Km/h. 

•  Bandwidths  of  1.5  to  20  MHz  depending 

•  Support  of  both  FDD  (Frequency  Division  Duplexing)  and  TDD  (Time  division  Duplexing) 

•  User  latency  of  5ms  and  control  latency  of  100  ms  going  from  idle  to  active 

•  Support  for  20  users  at  5MHz  bandwidth  and  400  users  at  higher  bandwidths. 

The  difference  between  the  LTE  and  prior  standards  is  the  reduction  of  latency  and  new 
modulation  scheme  that  has  improved  the  performance.  LTE  has  reduced  the  interfaces 
between  the  core  network  and  the  base  station  which  is  also  called  eNodeB.  In  LTE  eNodeB 
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does  the  both  jobs  of  mobility  management  and  radio  resource  management  which  reduce  the 
latency  between  these  tasks.  Figure  2-6  shows  the  structure  of  LTE  architecture. 

LTE 


c 


Core  network 


Figure  2-6  LTE  architecture 


LTE  is  divided  into  user  plane  protocols  and  control  protocols  which  all  reside  between  the 
eNodeB  and  the  UE  (user  end)  or  as  it  called  the  mobile  terminal  (also  User  Equipment).  These 
protocols  are  shown  in  Figure  2-7. 


User  Plane 


Control  Plane 


Radio  Resource 
Control  (RRC) 


Packet  Data  Convergence  Protocol 
(PDCP) 


Radio  Link  Control 
(RLC) 


Medium  Access  Control 
(MAC) 


Physical  Layer 
(PHY) 


Figure  2-7  User  plane  and  Control  plane  protocols 
The  protocol  layer  we  are  interested  in  is  the  Physical  layer  which  is  all  the  modulation;  coding 
and  antenna  schemes  are  located.  LTE  uses  OFDM  with  cyclic  prefix.  OFDM  is  a  type  of  multi¬ 
carrier  transmission  where  narrowband  channels  overlap  each  other  but  remain  orthogonal 
due  to  their  unique  orthogonal  basis  function.  Figure  3-3  shows  an  OFDM  carrier  division. 
Figure  2-8  shows  LTE  block  diagram. 
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Figure  2-8  OFDM  carrier  spectrum 
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Figure  2-9  LTE  Physical  layer  block  diagram 


The  two  important  blocks  that  we  are  interested  in  are  the  resource  mapping  and  the  antenna 
port  blocks.  Resource  mapping  has  the  task  to  assign  symbols  to  the  carriers  based  on  the 
quality  of  the  particular  channel.  LTE  transmission  is  based  on  radio  frames.  LTE  has  a  10  ms 
radio  frame  which  is  divided  into  10  sub-frame  and  each  1  ms  sub-frame  is  divide  into  2  slots  of 
0.5  ms.  Slot  is  the  smallest  unit  of  LTE  radio  structure  with  0.5  ms  time  duration  each  slot 
contains  6  or  7  symbols  for  normal  and  extended  cyclic  prefix  consequently.  The  length  of  cyclic 
prefix  is  related  to  the  channel  length  and  the  symbol  number.  A  larger  CP  is  intended  in  rural 
areas  and  the  shorter  CP  in  urban  areas  where  the  channel  length  is  smaller  and  higher 
throughput  without  a  loss  of  signal  orthogonality  is  required. 
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OFDM  symbols  have  either  normal  or  extended  cyclic  prefix.  In  case  of  a  normal  CP,  the  CP 
length  is  equal  to  5.2  us  for  the  first  symbol  and  4.7  us  for  the  remaining  symbols.  A  CP  length 
of  17.7  us  is  used  for  the  extended  CP  which  as  mentioned  before  is  used  in  areas  where  there 
is  a  long  delay  spread.  OFDM  symbol  duration  without  a  guard  interval  has  a  duration  of  66.7 
us.  Figure  2-10.  Shows  LTE  radio  frame  and  the  mentioned  timings. 


Radio  frame  (=  20  slots): 


Sub-frame  (=  2  slots): 


Slot  (=  6  or  7  OFDM  symbols): 


OFDM  symbol: 


0.5  ms 


* - > 

1  1  2 

3 

6  or  7 

<  > 

CP 

Figure  2-10.  LTE  radio  frame 


LTE  supports  channel  bandwidth  from  1.4  to  20  MHz  and  subsequently  subcarrier  sizes  from 
128  to  2048  carriers.  Sampling  time  of  it  has  been  carefully  placed  at  multiple  of  3.84  MHz  at 
30.72  MHz  for  back  compatibility  with  legacy  standard  like  UMTS  and  EDGE. 


The  subcarrier  frequency  space  is  divided  into  smaller  units  named  as  resource  block.  Each 
resource  block  consists  of  12  carriers  and  2  slots  which  can  virtually  be  shown  in  time  and 
frequency  domain.  A  resource  block  consists  of  84  resource  elements  for  normal  cyclic  prefix 
and  72  resource  elements  for  extended  cyclic  prefix.  A  resource  element  is  a  symbol 
transmitted  on  a  single  carrier  on  an  OFDM  symbol.  Carriers  in  a  resource  block  are  spaced  15 
KHz  apart,  therefore  an  RB  spans  the  bandwidth  of  180  KHz. 


In  order  to  obtain  channel  information  for  coherent  detection  and  synchronization  some 
carrier-symbol  cells  of  the  resource  block  is  dedicated  to  reference  signals.  These  reference 
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signals  are  either  Pilots,  which  are  known  symbols  transmitted  at  certain  RE's  or  dedicated 
channels  for  synchronization  such  as  PSS  (primary  sync,  series)  and  SSS  (secondary  sync,  series) 
channels.  Pilots  are  inserted  in  the  RB  depending  on  what  antenna  port  is  being  used.  In  a 
simple  single  antenna  scheme  pilots  is  placed  in 

One  downlink  slot  7^, 
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Figure  2-11  Resource  block  (RB)  and  Resource  Element  (RE) 


first  and  third  last  OFDM  symbol.  Pilots  are  placed  6  carriers  apart  in  frequency  domain  and  3 
symbols  apart  in  time  domain.  The  simplest  channel  estimator  scheme  needs  to  interpolate  in 
time  and  frequency  domain  to  obtain  the  channel  response  for  the  entire  resource  block. 
Obviously  more  sophisticated  methods  are  used  when  mobility  is  involved. 
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Figure  2-12  Resource  element  and  pilot  positioning  in  CP  for  two  antenna  ports 


In  a  multiple  antenna  scheme  pilots  are  positioned  such  that  they  are  orthogonal  to  other 
antenna  ports.  A  channel  estimation  algorithm  requires  calculating  the  channel  information  for 
each  of  the  antennas  separately.  In  a  multiple  antenna  scheme,  the  resource  elements 
occupied  for  the  pilot  in  one  antenna  port  is  empty  in  another  antenna  port  RE  to  avoid 
interference  with  each  other.  Figure  2-12  shows  the  pilot  positioning  for  a  normal  CP  scheme  in 
two  antenna  ports.  In  the  channel  estimation  section  we  will  show  how  the  pilots  are  used  for 
estimating  the  channel  information.  Before  we  proceed  into  detail  physical  layer  algorithms  we 
need  to  investigate  the  procedures  in  LTE  to  know  how  UE  and  base  station  communicate. 

Important  physical  layer  procedures 

The  most  important  procedure  that  UE  uses  to  initiate  a  communication  with  the  eNodeB  are: 


•  Cell  search 

•  Random  Access  procedures 


CELL  SEARCH 

Cell  search  is  the  method  that  UE  uses  to  find  a  suitable  base  station  to  communicate  with.  UE 
identifies  the  cell  which  is  associated  with  the  base  station  via  detecting  the  primary  and 
secondary  synchronization  signals  that  are  transmitted  from  the  base  station  on  regular  basis. 
When  discussing  the  synchronization  procedures  we  will  examine  the  cell  search  procedure  in 
more  detail. 
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What  important  for  us  is  that  cell  search  procedure  is  also  utilized  to  search  neighbouring  cells 
for  handover  procedures.  Resource  have  to  be  allocated  at  specific  time  to  minimize 
interference  with  cells  own  resources. 

RANDOM  ACCESS  PROCEDURES 

UE  requests  resources  from  base  station  using  the  RA  channels.  This  would  be  important  as  if 
utilized  extra  information  would  be  available  for  the  UE  to  estimate  its  position.  In  order  to 
initiate  such  requests  the  user  terminal  needs  to  send  known  preambles  which  are  the  Zadoff- 
Chu  sequences.  These  sequences  have  the  characteristic  that  they  build  an  orthogonal  basis 
each  and  when  transmitted  together  are  separable.  As  long  as  they  are  not  repeated  by 
different  users  there  would  be  no  collision.  When  base  station  detects  the  preamble  or  in  other 
words  the  RA  procedure  request,  the  terminal  can  send  further  information  and  initiate  the 
channel.  Random  access  channels  can  be  a  part  of  standalone  locationing  algorithm  as  it  utilizes 
the  eNodeB  and  UE  alone. 

In  order  to  conduct  the  experiments  as  realistic  to  the  LTE  physical  characteristics  we  briefly 
outline  the  LTE  characteristics. 

SUPPORTED  BW 

LTE  supports  the  following  bandwidths.  Each  bandwidth  then  would  have  specific  RB's  that  it 
can  support  which  in  turn  dictates  the  number  of  maximum  pilots  available  to  each  UE.  Table  2- 
1.  Shows  the  available  bandwidth. 


Channel  bandwidth  (MHz) 

1.4 

1.6 

3 

3.2 

5 

10 

15 

20 

Number  of  resource  blocks  with  FDD 

6 

N/A 

15 

N/A 

25 

50 

75 

100 

Number  of  resource  blocks  with  TDD 

6 

7 

15 

16 

25 

50 

75 

100 

Table  2-1.  Supported  Bandwidths 

LTE  DUPLEX 

LTE  transmits  in  FDD  or  TDD  full  and  half  duplex  method.  E-UTRA  TDD  band  40  is  the  band  that 
occupies  the  ISM  band  and  could  be  utilized  for  tests.  This  band  occupies  2300  to  2400  MHz. 
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LTE  SPECTRAL  MASK 

FCC  guidelines  for  the  LTE  transmission  masks  are  shown  in  Figure  8.  It  is  suggested  that  tighter 
requirements  be  placed  for  locationing  purposes  to  reduce  interference  from  neighbouring 
cells.  Possibly  this  should  go  with  better  pulse  shaped  OFDM  transmitter  and  receivers  for 
locationing  purposes.  This  is  feasible  if  non-commercial  radio  is  intended. 
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Figure  2-13  LTE  TX  spectral  mask 


OFDM  have  succeeded  as  a  high  data  rate  transmission  scheme  on  the  fact  that  it  combines 
high  data  rates  with  low  symbol  period  time.  Unlike  a  single  carrier  transmission  where  data 
rate  is  inversely  proportional  to  the  symbol  duration  Ts,  therefore  susceptible  to  channel  delay 
spread  Td,  OFDM  by  transferring  serial  data  symbols  in  parallel  increases  he  symbol  time  in 
expense  of  slight  frequency  bandwidth  increase. 

In  OFDM  stream  of  data  is  first  converted  from  serial  to  a  parallel  stream  and  then  transmitted 
onto  each  carrier.  The  transformation  of  serial  stream  into  N  parallel  data  streams  increases  the 
symbol  time,  which  is  the  time  it  takes  for  one  symbol  to  transmit,  by  a  factor  of  N.  Figure  9. 
shows  a  typical  OFDM  transmission  scheme. 
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High  symbol  rate 


Low  symbol 
rate 


Figure  2-14  OFDM  transmission 

As  seen  in  Figure  2-14  each  symbol  is  transmitted  via  one  carrier.  In  an  OFDM  system  total 
number  of  carriers  is  always  more  than  the  transmitted  symbols.  Some  of  these  extra  carriers 
are  used  as  guard  bands  and  some  used  to  transmit  reference  symbols.  Figure  3-10  shows  an 
OFDM  transmitter  and  receiver. 


Figure  2-16.  OFDM  transmitter  and  receiver 
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One  concept  in  modem  design  is  the  equalization  and  match  filtering.  OFDM  modulation 
scheme  has  solved  the  match  filter  concept  by  using  two  orthogonal  and  matched  basis 
function  which  is  sine  and  cosine  transformation.  An  iFFT  and  FFT  is  used  as  transmitter  and 
receiver  matched  filters.  In  such  scheme  equalization  of  the  signal  is  done  in  frequency  domain, 
therefore  prior  to  match  filtering  preserving  the  symbol  energy  in  one  FFT  block  is  very 
important. 

OFDM  solves  the  channel  effect  and  the  ISI  by  introducing  Cyclic  Prefix  and  by  simply  copying 
the  last  L  symbols  to  the  beginning  of  the  block  and  then  transmitting  through  the  channel.  By 
doing  so  any  channel  with  length  less  than  L  will  have  its  effect  preserved  in  the  CP.  In  this  way 
an  FFT  block  becomes  cyclo-stationary  and  it  would  be  unimportant  from  which  position  the 
FFT  block  is  taken  as  long  as  it  contains  the  whole  CP  as  a  part  of  the  block.  By  doing  this 
basically  the  effect  of  the  channel  is  reversed  and  OFDM  blocks  are  capable  of  equalizing  the 
channel  very  effectively  with  minimal  resources  provided  that  they  are  given  the  channel 
information  in  frequency  domain.  Figure  3-11.  shows  the  position  of  CP  in  an  FFT  block.  In  this 
figure  Tcp  is  the  duration  of  the  CP  which  depends  on  the  channel  it  is  going  to  be  used.  Tu  is  the 
actual  symbol  duration. 


Figure  2-17:  Cyclic  Prefix  in  FFT  block 

At  the  receiver  the  beginning  of  the  packet  is  found  and  from  that  position  and  FFT  block  is 
separated  and  an  ISI  free  block  is  passed  to  the  FFT  which  in  turn  transfers  back  the  signal  from 
what  is  known  time  domain  to  frequency  domain.  We  will  dig  into  the  math  of  such  operation 
next. 

Down  Link  Model 

Now  that  we  have  familiarized  ourselves  with  OFDM,  LTE  basic  facts  and  procedures  we  will 
discuss  the  LTE  downlink  specifically. 
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•  LTE  resources 


•  LTE  Synchronization 

•  LTE  Channel  estimation 

We  also  discuss  LTE  preambles  (Zadoff-Chu  sequence)  and  their  unique  characteristics. 

SYNCHRONIZATION 

No  resource  can  be  extracted  including  pilots  if  proper  synchronization  is  not  performed 
between  eNodeB  and  UE.  The  following  synchronization  has  to  be  performed: 

•  Carrier  frequency  synchronization 

•  Symbol  and  frame  timing  synchronization 

•  Sampling  clock  synchronization 

Without  the  above  three  synchronization  no  pilots  can  be  extracted  and  no  channel 
information  culd  be  calculated,  therefore  prior  to  channel  estimation  we  are  interested  in  how 
to  make  a  good  synchronization. 

CELL  SEARCH 

UE  (user  equipment)  needs  to  perform  a  cell  search  procedure  to  access  an  LTE  cell.  This 
involves  detecting  the  Primary  Synchronization  Signal  (PSS)  and  Secondary  Synchronization 
Signal  (SSS)  which  is  broadcasted  by  the  eNodeB  at  specific  radio  frame  locations  for  this 
purpose. 


Figure  2-18:  Cell  search  procedures 

As  seen  from  Figure  2-18,  in  order  to  examine  the  pilots  for  neighbouring  cell  and  associated 
cells,  SSS  and  PSS  has  to  be  performed  in  order.  In  a  radio  frame  the  PSS  and  SSS  signal  position 
is  different  for  FDD  or  TDD  transmission.  Figure  2-19  shows  the  position  of  these  signals  in  the 
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radio  frame  for  both  types  of  frames.  Prior  to  SSS  detection  the  CP  length  is  unknown  so  the 
synchronizer  blindly  should  search  for  the  SSS  in  the  given  locations.  Totally  4  different 
positions  has  to  be  searched  if  both  TDD  and  FDD  is  supported  in  an  algorithm  design.  PSS 
signal  always  proceeds  SSS  either  immediately  (FDD)  or  after  3  symbols  (TDD).  The  job  of  a 
locationing  algorithm  besides  extracting  information  from  the  channel  information  is  to  provide 
better  synchronization  algorithms  for  neighbouring  cells  if  it  intends  to  use  multiple  cell 
information  in  its  positioning  algorithms. 


We  mentioned  the  position  of  the  PSS  and  SSS  signals  in  time  domain  which  their  position  in 
the  radio  frames.  Now  we  are  interested  in  their  position  in  frequency  domain.  As  no 
bandwidth  information  is  known  prior  to  synchronization  the  PSS  and  SS  are  positioned  in 
centre  of  the  frequency  band  occupy  6  resource  blocks  (RB's)  which  is  the  minimum  number  to 
allocate.  The  PSS  and  SSS  have  62  symbols  which  are  all  mapped  to  62  central  carriers  in  the  6 
mentioned  RB's.  As  mentioned  before  each  RB  has  12  subcarriers  which  mean  from  72 
subcarriers  available  only  62  carriers  are  used  and  the  rest  is  not  used.  This  means  5  RE  at  each 
end  of  the  resource  block  is  not  used.  Figure  2-20  shows  the  position  of  the  RB's  and  the  empty 
RE's  at  the  lower  end  RB. 
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Figure  2-19:  (a)  FDD  frame,  (b)  TDD  frame 
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10  ms  Radio  Frame 


Figure  2-20:  PSS  and  RSS  position  in  RB's 

DETECTION  OF  PSS  AND  SSS 

PSS  uses  Zadoff-Chu  (ZC)  sequence.  ZC  sequences  are  Constant  Amplitude  Zero  Auto  correlation 
sequences.  Unlike  golden  codes  these  codes  are  non-binary  and  suitable  for  IQ  transmission. 
Equation  1  shows  the  odd  length  ZC  sequence: 

f-j27rg(n(n2+1))  +  t.n/  | 

/ Nzc\ 

aq(n)=e >  (1) 

Where  : 

qe{  1, ... ,  Nzc  —  1}  is  the  root  index. 

n  =  {0,1 . JVZC-1} 

ZeN  where  in  LTE  I  =  0  is  used. 

Rkk(  t)  =  Sn=0_1  ak(n)a*k(n  +  r)  =  S  (r)  (2) 

Where  R(.)  is  the  autocorrelation  function  of  ai<  at  position  r 
The  properties  of  these  sequences  are: 

•  ZC  sequences  have  constant  amplitude  in  time  and  frequency  domain. 
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•  ZC  codes  have  ideal  cyclic  autocorrelation  which  enables  the  user  to  detect  the  codes 
even  when  they  are  delayed.  If  the  code  is  designed  to  operate  at  delay  N.  The 
autocorrelation  of  a  delayed  sequence  will  generate  a  peak  at  the  delayed  time  as  long 
as  the  delay  is  shorter  than  the  designed  maximum  delay  of  N  samples. 

•  Multiple  orthogonal  sequences  can  be  generated  from  one  sequence.  The  cross 
correlation  of  two  sequences  have  a  constant  amplitude  of  1/^JNZC  if  and  only  if  the 
difference  of  sequence  indices  are  relatively  prime  to Nzc.  This  means  if  Nzc  is  chosen  as 
a  prime  number  Nzc  —  1  orthogonal  sequences  may  be  generated  by  just  changing  the 
indices  of  the  sequence. 

•  A  DFT  of  the  ZC  sequence  is  a  weighted  cyclic  shift  ZC.  This  means  ZC  can  be  generated 
in  frequency  domain  directly.  Later  some  papers  will  be  shown  that  use  this  property  to 
generate  ZC  sequences  efficiently. 


PSS  SIGNAL: 

As  mentioned  above  the  ZC  codes  can  be  generated  in  frequency  domain  directly.  The 
frequency  domain  codes  have  the  exact  same  characteristic  as  time  domain  ZC  codes.  A  PSS 
signal  uses  this  feature  and  transmits  a  63  length  ZC  code  on  63  carriers  with  middle  code 
punctured  to  avoid  dc  subcarrier.  LTE  uses  3  PSS  codes  for  3  possible  physical  layer  sections  in  a 
cell  which  means  one  PSS  code  for  each  group  in  the  cell. 

In  equation  (1)  q  is  the  root  of  the  ZC  sequence.  PSS  codes  are  defined  if  choosing  l  =  0  and  q  = 
29,  34,  25.  The  chosen  PSS  signals  have  low  sensitivity  on  frequency  offset  and  can  still  be 
detected  if  there  is  a  7.5  KHz  frequency  drift.  This  is  due  to  their  almost  flat  frequency  response. 
A  correlator  will  decide  which  code  has  been  transmitted.  The  received  signal  is  correlated 
against  the  existing  codes  and  the  position  of  the  peak  shows  the  offset  of  the  signal. 

SSS  SIGNAL: 

SSS  is  signal  is  an  M-sequence.  It  is  constructed  by  interleaving  two  31  length  BPSK  M-sequence 
codes  known  as  SSC1  And  SSC2.  SSC1  and  SSC2  are  cyclic  shifted  versions  of  a  31  length  M- 
sequence.  The  shift  of  each  code  is  mentioned  in  6.11.2.1-1  in  technical  specification  36.211  of 
3GPP  group.  SSC2  is  scrambled  by  a  sequence  depending  on  the  SSC1  shift  and  SSC1  is 
scrambled  based  on  the  PSS  signal  index. 
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SSS  detection  is  done  after  PSS  detection  and  therefore  the  system  has  the  channel  information 
at  this  point.  The  system  can  either  use  a  coherent  detector  or  a  non-coherent  one  which  is 
based  on  the  code  correlation. 

In  coherent  detection  the  performance  depends  on  the  presence  of  the  interferer  and  the 
quality  of  the  channel  estimation,  this  is  due  to  the  fact  that  the  interfere  uses  the  same  PSS 
signal  therefore  the  SCI  and  SC2  signals  would  be  identical.  In  a  non  coherent  detection  the 
performance  depends  on  the  coherent  bandwidth  of  the  channel  which  should  be  greater  than 
the  6  resource  blocks  being  used  by  the  SSS  signal.  This  ensures  that  the  SS  signal  would  not  be 
changing  during  the  transmission. 


Channel  estimation 

A  locationing  algorithm  needs  channel  estimation.  A  coherent  detector  is  used  that  can  utilize 
the  phase  and  amplitude  of  the  received  signal.  A  common  way  to  estimate  the  channel 
information  is  to  use  a  known  signal  such  as  a  pilot.  These  known  signals  are  called  Reference 
Signals  (RS's).  One  can  realize  an  OFDMA  signal  as  a  three  dimensional  frame  of  time,  frequency 
and  space  signal.  RS  signals  are  placed  in  the  LTE  frame  at  known  time  and  frequency  positions 
on  an  optional  or  chosen  space  which  is  the  antenna  port.  The  position  in  frequency  is  called 
subcarrier  and  the  position  in  time  is  known  as  slots  and  the  smallest  element  in  this  two 
dimensional  domain  is  called  the  Resource  Element  (RE). 

LTE  downlink  uses  5  different  reference  signals  (RS)  based  on  the  LTE  release  number. 

•  Common  reference  signal  which  are  available  to  all  UEs  (Pilots) 

•  Demodulation  reference  signal  specific  to  the  UE 

•  Multimedia  broadcast  Single  Frequency  Network  MBSFN  specific  reference  signal 

•  Positioning  reference  signal  (Release  9  onwards) 

•  Channel  state  Information  reference  signal  (Release  10  LTE  advanced) 

Reference  signals  are  transmitted  from  eNodeB.  LTE  introduces  antenna  ports  as  ports  to 
transmit  RS  which  can  either  be  physical  antenna  or  virtual  ports.  From  UE  point  of  view  any 
measurement  obtained  from  these  antenna  ports  defines  the  characteristic  of  that  port. 
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What  we  are  going  to  first  target  for  the  channel  estimation  is  the  common  reference  signals  or 
the  cell  specific  signals.  These  RS  signals  are  available  to  all  UE's  in  the  cell  and  it  is  a  standard 
RS  since  release  8,  these  are  the  primary  RS  for  channel  estimation  purposes. 

In  Figure  15  an  R0  indicates  a  RS  for  antenna  port  zero.  The  transmitted  RS  signal  is  a  QPSK 
signal  defined  in  equation  (3). 

n,ns=  ^[1  -  2c(2m)]  +;^[1  -  2c(2m+  1)]  (3) 

In  equation  (3),  m  is  the  reference  signal  index,  ns  is  the  slot  number  and  l  is  the  symbol 
number.  C(m)  is  taken  from  a  length  31  gold  sequence  which  is  initialized  based  on  the  RS  type. 
The  initialization  is  done  every  OFDM  symbol  and  the  value  depends  on  the  cell  identity,  Nfp11. 
Each  cell  uses  a  cell  specific  frequency  shift  equal  to  Nfp11  mod  6  which  avoids  frequency 
collisions  between  the  RS  signals  from  neighboring  6  cells.  This  is  an  important  feature  when 
using  the  channel  estimation  from  neighboring  cells  to  calculate  the  location. 
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Figure  2-21:  Cell  specific  RS  for  normal  CP  length 
We  would  leave  the  rest  of  the  RS  signals  for  now  till  later  phases  of  the  project  but  before 
leaving  this  section  we  need  to  mention  the  next  important  RS  which  is  the  UE  specific  RS.  UE- 
specific  RS  are  used  only  between  eNodeB  and  particular  UE  and  are  transmitted  as  of  any 
other  PDSCH  (Physical  Downlink  Shared  Channel)  signal.  The  difference  is  that  these  RS  signals 
are  mapped  to  specific  RB's  that  target  the  specific  UE  and  therefore  contain  UE  specific  RS 
signals.  This  can  be  beam-forming  signals  or  extra  RS  signals.  We  will  go  into  details  of  these 
specific  RS  signals  later. 

We  will  consider  the  channel  h(r,  t)  to  be  the  channel  time  varying  complex  signal  at  delay  r 
and  time  instant  of  t.  In  a  discrete  form  this  channel  is  modelled  as  h[l,  kT]  with  l  as  the  delay 
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at  time  sampling  instant  kT.  We  also  consider  the  channel  to  be  approximated  by  an  FIR  vector 
of: 

h[k]  =  h[0,k] . h[L-l,k]  (4) 

where  L  is  the  channel  delay. 

Practical  channel  estimation  for  LTE  is  done  in  frequency  domain.  If  the  received  signal  is: 
y  =  h  *  x  +  n  (5) 

The  FFT  of  the  received  signal  is: 

Y  =  HX  +  N  (6) 

Where  X  is  the  Nxl  FFT  of  the  transmitted  symbols  and  FI  is  the  LxN  cyclic  FFT  of  the  channel 
information  and  N  is  a  zero  mean  complex  cyclic  white  noise  with  covariance  of  C.  In  order  to 
obtain  the  channel  matrix  FI  it  is  sufficient  to  correlate  the  X  with  the  known  transmitted 
symbols  and  obtain  FI: 

Y  =  H(X.Xh)  +  NXh  =  H  +  z  =  H  (7) 

If  X  is  the  reference  symbols  and  has  zero  mean  and  unit  energy  then  z  is  still  a  circular  zero 
mean  noise  and  the  result  would  be  the  NxL  estimated  channel  where  L  is  the  channel  length. 

As  shown  in  Figure  15,  the  channel  estimation  would  result  in  estimating  the  channel  for 
locations  in  frequency-time  lattice  that  the  RS  signal  is  present.  In  order  to  estimate  the  channel 
for  the  rest  of  the  subcarriers  approximate  estimation  of  the  channel  at  the  frequency  and  time 
lattice  is  required.  If  the  frequency  and  time  is  considered  statistically  independent  then  we  can 
divide  the  estimation  in  2  domains  as  a  single  domain  linear  estimation  which  in  its  simplest 
form  will  be  interpolation. 

More  sophisticated  approaches  to  obtain  the  two  dimensional  estimates  would  be  LS  least 
square  and  MMSE  approach.  Any  approach  in  time  domain  has  been  proven  to  be  impractical 
although  from  estimation  point  of  view  it  would  have  resulted  in  more  accurate  unbiased 
estimations.  Chapter  8.4  in  [1]  compares  different  channel  estimation  methods  but  the  most 
practical  channel  estimation  for  an  immobile  user  is  the  interpolation  method  described  above. 
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Simulation  Environment 


The  system  requires  TOA  and  channel  estimation  from  multiple  eNodeB's  at  the  UE  position.  As 
mentioned  in  previous  sections  the  synchronization  methods  provided  in  LTE  system  would 
provide  the  TOA  of  signals  from  different  cells  to  the  UE  and  the  channel  estimation  methods 
described  in  previous  section  would  provide  the  channel  information  needed.  In  order  to  set  up 
such  simulation  environment  a  single  antenna  downlink  scheme  as  shown  in  Figure  16  is 
provided: 


code  words  layers  antenna  ports 


Figure  2-22.  Downlink  Physical  Chain  per  eNodeB 


The  OFDM  signal  from  different  eNodeB  goes  through  a  Lagrange  interpolator  to  change  phase 
according  to  their  distance  to  the  UE  and  then  multiplexed  with  other  eNodeB's  and  passed 
through  the  channel.  Figure  2-23  shows  an  alternative  that  the  signal  is  individually  passed 
through  the  channel  and  delay  is  performed  on  the  signal  after  the  channel. 


Up  sampled  K  > 


Course  chip 
delay 


>  Lagrange  interpolator  fine  delay 


Figure  2-23.  Downlink  channel 


On  the  receiver  the  PSS  and  the  SSS  signal  is  extracted  and  the  delay  associate  with  the 
intended  eNodeB  is  considered.  The  interpolated  signal  for  each  intended  eNodeB  is  passed  to 
their  channel  estimator  and  the  channel  is  extracted  for  each  eNodeB-UE  path. 
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The  LTE  observables  which  are  relevant  are  essentially  the  known  synchronization  signal 
segments  that  are  part  of  the  resource  blocks.  While  it  is  a  reasonable  assumption  that  the  LTE 
AN  signals  are  frequency  synchronized  (which  can  be  accurately  obtained  from  GPS  signals 
amongst  other  potential  sources),  it  is  not  a  reasonable  assumption  in  many  cases  to  assume 
that  the  AN's  are  accurately  time  synchronized.  This  depends  as  the  network  may  be 
specifically  calibrated  for  SGL  purposes.  Tight  timing  synchronization  amongst  the  AN's  is  not 
necessary  for  optimum  data  communication  functions  nor  is  it  necessary  for  facilitating 
network  based  location,  hence  there  is  little  impetus  for  the  service  provided  to  accurately  time 
synchronize  the  network.  Consequently,  the  SGL  processing  developed  will  consider  the  two 
modes  when  the  LTE  network  is  tightly  time  synchronized  and  when  it  is  not. 
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Chapter  3  CV  observables  and  processing 


3.1  Overview 

The  SGL  observables  created  from  LTE  range  measurements  are  useful  as  the  wireless  signals 
are  generally  ubiquitous  throughout  the  environment  that  the  UE  is  immersed  in.  The  LTE  AN's 
are  generally  frequency  synchronized  but  may  not  be  time  synchronized  and  with  possible 
uncertainty  in  terms  of  the  AN  transmitter  location.  Furthermore  the  LTE  signals  are  subject  to 
a  significant  amount  of  multipath  which  for  isolated  SGL  estimates  renders  them  virtually 
useless.  However  the  SLAM  algorithm,  developed  in  the  previous  chapter  estimates  the 
trajectory  of  the  UE  in  light  of  all  of  these  uncertainties  partially  ameliorating  the  detrimental 
uncertainty  due  to  multipath.  CV  observables,  complement  the  LTE  observables  in  that  they 
provide  high  accuracy  over  shorter  segments  of  the  trajectory.  Hence,  relative  updates  in  the 
UE  position  can  be  achieved  based  on  CV  observables.  They  also  have  the  very  important 
attribute  of  being  independent  of  any  LTE  signalling  providing  observables  where  LTE  downlink 
signals  are  not  available  or  compromised  in  some  way.  Additional  notable  advantages  of 
incorporating  CV  is  that  it  is  possible  to  detect  if  the  LTE  signals  are  LOS  or  NLOS  over  small 
trajectory  segments.  Detection  of  a  small  LOS  trajectory  segment  implies  that  beamforming 
can  be  used  on  the  LTE  signals  which  are  highly  effective  in  determining  the  relative  position  of 
the  UE  with  respect  to  the  AN.  Another  advantage  of  using  CV  is  that  these  observables  are 
immune  to  LTE  signal  jamming  or  spoofing.  Finally  CV  observables  can  be  sourced  from  small 
uncalibrated  cameras  which  are  typically  part  of  the  UE  handset  such  that  additional  hardware 
to  generate  the  CV  observables  is  negligible. 

The  basic  method  of  generating  CV  sensor  inputs  is  to  track  features  on  the  ground  manifold 
surface  updating  outputs  of  the  differentials  in  the  pose  and  the  translational  motion  relative  to 
the  ground  manifold  at  the  frame  rate  interval  of  the  camera  as  illustrated  in  Figure  3.1.  The  CV 
update  rate  of  between  10  to  30  frames  per  second  (fps)  is  possible  while  maintaining  real  time 
processing  which  is  sufficiently  fast  to  provide  highly  accurate  relative  updates  of  the  UE 
position  and  orientation.  Ego-motion  is  essentially  the  process  of  determining  the  relative 
position  and  orientation  of  the  camera  based  on  what  it  observes.  What  is  observed  is  the 
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surrounding  3D  scenery  in  the  world  reference  that  is  mapped  onto  the  image  plane  of  the 
camera.  If  features  of  known  identity  and  location  (relative  to  the  world  reference)  can  be 
observed  then  we  have  ego-motion  in  the  absolute  sense.  However,  a  more  common  scenario 
is  that  the  features  cannot  be  uniquely  classified  or  do  not  have  a  location  that  is  known 
absolutely.  Hence  only  relative  motion  and  orientation  of  the  camera  is  possible.  In  the  latter 
case,  it  is  generally  assumed  that  the  feature  points  are  static  relative  to  the  world  reference  and 
hence  the  relative  ego-motion  is  possible.  In  using  these  static  feature  points  of  unknown 
absolute  location,  the  ego-motion  algorithm  actually  simultaneously  determines  its  own 
trajectory  relative  to  the  array  of  feature  points  such  that  all  the  observations  are  consistent. 

That  is,  in  the  end,  the  relative  positioning  of  both  the  trajectory  and  fp’s  is  then  established 
which  is  a  form  of  SLAM,  (simultaneous  location  and  mapping).  If  a  set  of  known  fp’s  with 
known  locations  in  the  world  reference  is  observed  by  the  camera  then  the  combined  array  of 
fp’s  and  trajectory  can  then  fixed  absolutely  to  the  world  reference.  Hence  absolute  mapping 
and  location  is  achieved. 

Feature  points  that  are  suitable  for  tracking  are  determined  by  the  GF2T  (Good  Features  to 
Track)  and  then  tracking  them  frame  to  frame  with  LKP  (Lucas  Kanade  Pyramid)  method.  The 
LKP  is  used  to  establish  the  necessary  correspondence  between  the  fp’s  in  different  frames 
such  that  the  motion  of  the  fp’s  can  be  determined.  The  set  of  fp’s  used  will  be  assumed  to  be 
stationary  relative  to  the  world  reference  frame.  We  can  assume  different  sets  of  constraints. 

For  instance,  we  can  assume  that  the  FOV  of  the  camera  is  confined  to  the  floor  surface  or  the 
ceiling  of  a  building.  We  may  assume  that  the  fp’s  are  distributed  instead  in  the  3D  space  but 
perhaps  confined  to  a  set  of  planes.  These  planes  may  be  known  to  intersect  at  90  degrees  or 
may  be  more  random.  Generally  we  will  find  that  the  camera  provides  a  wealth  of  observations, 
however,  the  ego-motion  problem  is  generally  ill-conditioned.  Hence  while  it  is  possible  to  do 
full  6  DOF  ego-motion  (full  evaluation  of  T  and  R)  with  a  single  camera,  two  cameras  or  a 
projector  and  camera  is  in  the  Kinect  application  is  typically  required  to  make  the 
implementation  robust. 

We  will  start  the  ego-motion  development  with  a  highly  constrained  motion  in  2D.  The  2D 
problem  will  involve  tracking  of  stationary  fp’s  coincident  with  a  2D  manifold  surface  such  as  the 
floor  surface  or  a  wall.  As  is  illustrated  in  Figure  3-1 ,  the  camera  is  oriented  such  that  it  is 
looking  at  a  flat  2D  manifold  surface  which  could  be  a  ground  or  floor  surface.  As  the  camera  is 
moved  across  this  surface,  there  will  be  a  flow  of  features  in  the  field  of  view  which  is  related  to 
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the  motion  of  the  camera  itself.  This  2D  problem  is  significantly  simpler  than  the  3D  case  where 
the  camera  has  to  perceive  a  3D  environment,  usually  with  multiple  cameras  or  correlating 
multiple  views  of  the  same  camera.  Throughout  this  chapter  we  will  assume  a  perfectly 
calibrated  camera  that  behaves  as  an  ideal  pin-hole  camera  which  is  reasonable  as  camera 
calibration  can  be  assumed  as  it  is  a  robust  procedure. 


Figure  3-1  Ego-motion  relative  to  a  2D  manifold  surface 


3-2  System  Definitions 


In  this  section  a  set  of  necessary  definitions  will  be  given  that  will  facilitate  the  development  of 
the  ego-motion  algorithms.  Initially  only  a  single  camera  will  be  considered.  In  each  image 
frame  made  available  from  this  camera  as  it  moves,  there  is  a  set  of  observable  fp’s.  New  fp’s 
are  typically  identified  by  a  corner  feature  detector.  We  will  denote  this  as  the  Good  Features 
two  Track  (GF2T)  algorithm.  The  features  are  matched  from  one  frame  to  the  next  to  determine 
the  correspondence  based  on  the  Lucas  Kanate  Pyramid  (LKP)  routines. 

The  set  of  fp’s  generated  from  these  algorithms  for  the  tth  frame  is  denoted  as  {/**,//,,}  ■  Here 

the  index  t  denotes  the  frame,  k  denotes  the  feature  index  and  {x,y}  denote  the  integer  indices 

of  the  location  of  the  feature  point  as  referenced  to  the  image  plane  of  the  camera.  The 
analysis  is  primarily  focussed  on  a  monocular  ego-motion  single  camera  video  frame  sequence 
consistent  with  what  is  in  the  handheld  device.  We  will  start  with  Ego-motion  for  a  2D  plane  as 
shown  in  Figure  3-1  followed  by  ego-motion  in  a  3D  environment.  Differential  ego-motion  will  be 
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considered  which  is  based  on  determined  by  how  the  fp’s  have  moved  from  the  previous 
camera  image  to  the  current  image  as  captured  in  the  fp  data  and  v ,  fky  j 

respectively.  From  this  a  differential  or  incremental  update  can  be  determined  by  relating  this  to 
the  underlying  transformation.  As  established  in  the  previous  chapter  on  optical  flow,  finding  the 

correspondence  between  |/^,/^}  and  }  is  not  trivial.  Undetected  errors  made  in 

the  fp  correspondences  will  result  in  significant  distortions  of  the  estimated  differential  motion 
update.  Parameters  for  the  LKP  routine  have  to  be  carefully  selected  to  obtain  robust  fp 
correspondence  assignment.  In  the  following  chapter  on  Bayesian  filtering  emphasis  will  be 
given  to  the  overall  uncertainty  of  the  fp  migration. 

The  ego-motion  geometry  will  be  based  on  the  world  centered  and  camera  centered  coordinate 
systems  as  described  in  conjunction  with  the  perspective  transformation  in  a  previous  chapter 
with  the  following  definitions: 

Oc  -  camera  origin  in  world  coordinate  space 

Ow  -  world  origin  in  world  coordinate  space  (Usual  definition  isOw  =[0,0,0]r  ) 
flk  -  fp  in  the  world  coordinate  frame  with  coordinates  of  {fwXx,flXy,flXz} 
f‘k  -  fp  in  the  world  coordinate  frame  with  coordinates  of  {flk>x,flk>y,flk>z\ 

{XW,YW,ZW}  -  world  coordinate  system  directional  unit  vectors  (left  hand  system) 

{XC,YC,ZC}  -  camera  coordinate  system  directional  unit  vectors  (left  hand  system) 
n  -  image  plane,  given  by  the  plane  Z c  =  / 

{x,y}-  unit  vectors  of  2D  camera  image  plane  n 

Figure  3-  2  is  an  illustration  of  the  mapping  of  a  fp  between  the  world  and  camera  coordinate 
systems.  It  is  important  to  realize  that  the  rays  between  fc  k  and  f{v  k ,  as  illustrated  in  Figure  3- 
2  are  essentially  lines  of  ambiguity  where  the  fp  actually  originates  from. 
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camera 


Figure  3-2  Projection  of  fp  in  the  world  coordinate  system  onto  camera  2D  image  plane 
The  projective  transformation  defined  earlier  resulted  in  the  transformation  of 


P,=R(P„-T)  (1) 

relating  the  position  vector  in  the  world  coordinate  frame  to  that  of  the  position  vector  of  the 
camera  coordinate  frame.  Here  T  is  the  translation  vector  T  =  Oc  -Ow  and  the  rotation  matrix 

is  given  as  the  projection  coefficients  of  the  unit  vectors  of  {XW,YW,ZW}  onto  {Xc,Yc,Zc}as 


R  = 


xwnxc 

xwcyc 

x_tz. 


Ywnxc  zwnxc 

YW®C  Zw^c 

ywdzc  zwnzc 


(2) 


where  u  denotes  the  dot  product  between  the  unit  vectors.  (1 )  relates  the  fp’s  referenced  to  the 
world  and  camera  references  as 


(3) 


Expressing  R  in  terms  of  row  vectors 
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such  that  the  homogeneous  camera  coordinates  are  then 
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From  which  the  non-homogeneous  image  plane  coordinates  are  written  as 

,R,(p.-T) 

x  =  f - 7 - f 

R,(P.-T) 

_  .R,(P.-T) 

J  ^(p.-t) 

Z  =  f 


(6) 


In  the  application  of  current  interest,  a  handheld  device  is  hovering  over  a  2D  flat  ground 
surface  such  that  the  fp’s  are  residing  on  the  surface  the  surface  of  the  world  reference  as 

Zw  =  0  as  in  figure  1  such  that  fw  k  z  =  0 .  We  also  assume  that  the  fp’s  are  stationary  such  that 
for  the  world  coordinates  we  can  drop  the  “t”  index.  Then  we  have 


f 

J  c,k,x 

fw,k,x 

/•' 

*7  c,k,y 

=  H 

fw,k,y 

1 

where  H  was  defined  earlier  as 


X  hu  hx  ; 

tel  tel 

H  = 

7^2.i  7/22  H23 

= 

[RJ,  tel  -Rt 

h31  h32  h3  3_ 

tel  tel 

(7) 


(8) 


The  nonhomogeneous  coordinates  are  then 


(9) 


f*  =  f*  /  ff 

J  k,x  J  c,k,x  /  J  c,k,z 

f  =  f  /  f* 

J  k,x  J  c,k,y  /  J  c,k,z 


Hnfw,k,x+H32fw,k,y+H33 

H,,UX  +  H \lfw,k,y  +  ^33 


In  many  cases  of  practical  interest  where  the  fp’s  as  distributed  such  that  fc  k  z  does  not  vary 
significantly  then  the  weak  perspective  assumption  can  be  used  such  that  fc  k  z  is  taken  to  be  a 
constant  denoted  asc0  .  Then  we  have 


cJL=Hnf:Xx+Hl2f^y+Hl3 

(10) 

c0Kx=H2j:Ax+H22f:Ay+H23 

Note  that  we  are  really  implying  that  f‘kz  is  constant  relative  to  both  time  and  the  fp  index  k. 

With  this  the  perspective  matrix  reduces  to  that  of  the  affine  transformation  which  is 
conveniently  redefined  as 


f 

a  b  c 

fw  ,x,k 

J  x?/c 

- 1 

_ i 

fw,y,k 

1 

where  { a,b,c,d,e,f }  are  coefficients  dependent  on  the  geometry  of  the  camera  orientation  and 

position  relative  to  the  world  reference.  This  weak  perspective  transformation  or  affine 
transformation  is  clearly  applicable  where  the  camera  axis  Zc  is  collinear  with  the  Zw  axis. 

That  is  the  camera  can  translate  and  rotate  in  azimuth  but  have  not  tilt.  The  affine 
transformation  is  still  a  reasonable  approximation  for  cases  where  there  is  a  slight  tilt  provided 
that  the  fp’s  considered  in  the  world  plane  of  zw  =  0  are  close  to  the  intercept  point  of  the  Zc 

axis  in  the  zw  =0  plane. 

Finally  define  At  as  the  transformation  corresponding  to  the  tth  frame  such  that 
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(12) 


ft 

fw,k,x 

J  k,x 

=  At 

fw,k,y 

J  k,y 

1 

3.3  2D  translation  ego-motion 

Consider  the  camera  oriented  horizontally  over  the  Zw  =  0  plane  as  illustrated  in  figure  3  such 
that  the  perspective  rotation  matrix  is  R  =  I3 .  Denote  the  perspective  translation  vector  as  T, 
where  the  subscript  t  denotes  the  frame.  For  the  2D  ego-motion  case  at  a  constant  height  h  we 
have  T,  =  y‘T  -hj  where  j  is  the  displacement  along  the  x  and  y  axis  of  the 

camera  frame  relative  to  the  world  reference  and  h  denotes  the  height. 

Begin  with  the  simplified  1 D  case  of  motion  only  along  the  Xw  axis  such  that  Tt  =  0  -/zj 

Assume  that  the  sets  of  feature  points  of  f‘ck  for  the  current  frame  and  f'~k  from  the 
previous  frame  are  available. 


Figure  3  Simplified  ID  translation  in  the  Xw  direction 
Based  on  this  definition  the  affine  transformation  is  given  as 
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(13) 


1  0  -Xj 

0  1  0 

Typically  the  scaling  of  f  jh  is  ill  defined  as  the  focal  length  of  the  camera  is  typically  not  well 

known.  Hence  we  replace  this  with  a  scaling  factor  of  s  which  is  assumed  to  be  known  from  the 
camera  calibration.  The  overlay  of  the  kth  fp  in  two  successive  frames  is  illustrated  in  Figure  4. 

Note  that  we  can  determine  the  differential  change  AAt  □  At  -  At_,  but  we  cannot  determine  the 
absolute  At .  From  (12)  we  have 


ft 

ft- 1 

J  k,x 

J  k,x 

ft 

ft "I 

J  k,y 

J  k,y 

fw,k,x 

fw,k,x 

(A,-Aw) 

fw  ,k,y 

=  AAt 

fw,k,y 

1 

1 

(14) 


Based  on  the  case  where  we  have  motion  only  along  the  x  axis 


AAt 


=  s 


0  0  —Ax‘t 
0  0  0 


(15) 


where  Ax‘T  □  x‘r  -x‘T  1  and  s  is  a  scaling  factor.  Based  on  similar  triangles  we  can  derive 

~f£)  (16) 

S  V  7 
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fw,k,x 

Figure  4  Similar  triangles  used  to  determine  the  geometry  between  the  camera  image  plane 
and  ground  plane 


Note  that  we  can  estimate  the  relative  1 D  ego-motion  based  on  a  single  fp.  If  the  fp  is  a 
landmark  or  anchor  point  then  we  can  determine  the  absolute  ego-motion  from  just  a  single 
frame  with  the  assumption  that  the  camera  is  precisely  calibrated.  As  stated  earlier,  the 
practicality  of  CV  ego-motion  estimation  is  that  the  relative  camera  motion  can  be  determined 
from  arbitrary  and  previously  unknown  fp’s  that  are  only  assumed  to  be  temporally  static  relative 
to  the  world  (ground  manifold)  reference  frame.  An  issue  with  fp’s  is  that  they  are  subject  to 
random  fluctuations  from  one  frame  to  the  next  as  observed  in  the  previous  chapter.  Typically 
the  random  fluctuation  is  independent  from  one  fp  to  the  next.  Hence  the  random  deviation  of 

Ax‘t  can  be  reduced  by  using  the  complete  set  of  observed  fp’s  in  the  cameras  FOV.  This 
results  in  an  over-determined  set  of  constraints  for  Ax,  which,  in  the  absence  of  further  statistical 
information,  can  be  optimally  used  in  a  least  squares  sense  as 


fu-fi 


t- 1 


t- 1 


flC,x  flC,x 


j^Ax^ 


(17) 
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The  fp  constraint  relations  are  expressed  in  this  rather  cumbersome  way  to  establish  the 
notation  that  we  will  require  shortly  when  we  have  more  than  one  ego-motion  variable.  Note 
that  we  have  sets  of  over-determined  equations  of  the  form  of  M  =  AP  where  M  denotes  the 
measurement  vector,  A  the  model  matrix  and  P  the  parameter  vector  of  the  unknowns.  The 

least  squares  solution  is  P  =  ( A7A)r  ArM .  This  approach  will  be  used  for  all  of  the  ego- 

motion  calculations.  In  this  case,  Ax  becomes  the  average  of  the  data  contained  in 


M  = 


(18) 


ft  _  ft- 1 
_JK,x  J K,x  _ 

Next  consider  the  more  general  motion  of  a  translation  in  both  the  xw  and  yw  directions  from 
the  frame  t-1  to  the  current  frame  t .  Use  the  notation  that  we  have  developed.  We  now  have 
T;  =  [4  y\  -hj  such  that  the  affine  transformation  is  given  as 


1  0 


0  1 


(19) 


Again  assume  that  the  sets  of  feature  points  of  fc  k  for  the  current  frame  and  f‘k[  from  the 
previous  frame  are  available.  The  differential  change  AAt  □  At  -  At_j  is  now 


AA 


t 


0  0  -Ax't 
0  0  -A y‘r 


(20) 


where  Ax‘T  □  x‘r -x‘T  1  and  A y‘T  □  y‘T -y‘T  1  .  Extending  the  notation  developed  for  the  ID  case 


ft  _  ft- 1 

Jk,  x  J  k,x 

'1  O' 

sAx‘t 

ft  _  ft-\ 

J  k,y  J  k,y 

_0  1_ 

.s’ Ay  j 

This  is  of  course  separable  into  two  independent  sets  of  constraints  over  the  K  fp’s  as 
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(22) 


These  are  solvable  in  the  LS  sense  which  in  this  case  just  amounts  to  determining  the  mean 
(over  k)  as 


-sAxt=meaniflx  -/£) 

a  (,t  ,m\  (23) 

-sAyt=mean[fky-fky) 

This  provides  the  estimate  of  the  differential  motion  of  the  camera  based  on  two  consecutive 
frames.  If  the  requirement  was  essentially  to  only  estimate  the  incremental  translation  of  the 
camera  between  the  two  frames  with  no  prior  information  then  (17)  would  constitute  the  best 
estimate  possible.  However,  we  are  generally  interested  in  the  translation  motion  over  multiple 
frames.  Extending  the  LS  solution  over  multiple  frames  significantly  improves  the  performance. 
For  example  if  we  make  the  additional  assumption  that  the  translation  velocity  is  constant  over 
three  frames  then  we  can  extend  the  LS  solution  as 


rt  _  ft-l 
J  K,x  JK,x 

rt+l  rt 

J l,x  ~  J l,x 

rt+l  rt 

J K,x  ~  J K,x 


fly 


7-1 


'fly 


ft  _  ft-\ 

J  K,y  J  K,y 
rt+l  rt 

Jl,y  ~  Jl,y 

rt+l  rt 

J K,y  ~  J K,y 


(24) 


Note  that  generally  this  would  result  in  a  variance  in  the  estimate  of  Axt  and  A yt  that  is  half  of 
the  variance  of  using  the  single  pair  of  frames.  The  message  here  is  that  if  the  camera 
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trajectory  is  smooth  then  a  significant  reduction  of  the  variance  of  the  ego-motion  variable 
estimation  is  possible. 

An  example  of  1  Dmotion  of  a  translating  square  subjected  to  noise.  The  first  frame  of  the 
video  is  shown  in  Figure  5  together  with  the  four  tracked  fp’s  as  green  donuts  and  the  LKP 
trajectory  of  each  fp  shown  in  blue.  White  uniform  noise  is  then  added  to  each  frame  that  is 
independent  from  frame  to  frame.  The  image  is  then  smoothed  which  gives  the  noise  a  blotchy 
texture. 

1 

Figure  3-5  Tracked  square  with  fp’s  indicated  and  the  LKP  trajectories 

The  plot  in  Figure  3-6  is  the  tracked  fp’s  which  is  very  clean  as  the  motion  in  x  is  3  pixels/frame 

which  is  an  integer.  Figure  3-7  shows  the  differentials  of  Ax,  and  A y, .  The  velocity  in  x  is 

changed  to  2.5  pixels/frame.  Figure  3-8  shows  the  differential  of  Ax,  which  now  has  a  ripple  due 
to  the  video  quantizing  effect.  To  avoid  this  problem  it  is  necessary  to  interpolate  the  video. 
Figure  9  shows  Ax,  and  Ay,  for  vx=2.5  and  vy=1 .2  showing  the  quantization  effects  of  the  block 
pixel  movement. 


X 


Figure  6  Trajectories  of  fp’s 
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Figure  3-7  Recovered  differentials  of  Ax,  and  Ay,  with  vx=3 


Figure  3-8  Recovered  differentials  of  Ax,  and  A yt  with  vx=2.5 
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Figure  3-9  Recovered  differentials  of  Ax,  and  Ay,  with  vx=2.5  and  vy=1.2 
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3-4.  Rotation  of  the  World  axis 


Next  consider  the  camera  at  a  fixed  position  such  that  x‘T  =  y‘T  =  0 .  As  before  the  height  is  fixed 
at  z‘T  =-h.  The  camera  is  rotated  in  a  positive  sense  around  the  zw  axis  of  (oz)  as  illustrated 
in  Figure  10  such  that  the  affine  transformation  is 


xcnxw  xenrw 

YcnXw  Ye[Yw 


0 

0 


cos(az) 

sin(az) 


-sin(az)  0 
cos(az)  0 


(25) 


Note  that  we  define  the  sense  of  rotation  relative  to  the  right  hand  but  we  assume  a  left  hand 
coordinate  system  for  both  the  camera  and  the  world  as  before. 


Figure  3-10  Rotation  offp  in  camera  FOV  by  az 
As  before,  we  have 
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=  At 
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1 

As  there  is  no  translation  this  transformation  can  be  simplified  to 
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(27) 


where  a[  is  the  azimuth  rotation  angle  between  the  world  and  camera  frame  at  time  t  which  we 
cannot  determine  absolutely.  However,  we  can  determine  relative  changes  as 
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where 

A  a‘z  =  a\  —  a'f 

For  notational  convenience  write 


C,  =  cos^Aa' ) 
St  =  sin^Aa' ) 

such  that 


(29) 


K  = 


C,  -s, 
st  a 


(30) 


73/143 


and  then  fk  =  R specifically  for  the  kth  fp.  Note  that  for  pure  differential  rotation  at  a 
constant  height  the  scaling  factor  cancels  out.  Combining  all  of  the  K  fp’s  in  a  LS  formulation  to 
estimate  Aa‘z ,  it  is  most  convenient  to  assume  the  redundant  linear  estimation  of 
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The  total  accumulated  rotation  can  be  inferred  from  the  rotation  matrix  of 


(31) 


/;=riR£//  P2) 

p=\ 

Hence  by  multiplying  the  estimated  differential  rotation  matrices  we  can  determine  the  total 
azimuth  rotation  angle  from  the  initial  frame  to  the  tth  frame.  This  is  approximately  equivalent  to 
the  processing  of 


A az  t  =  atan 2  ( Ct ,  Sf ) 

t 

a_=}  A  a,  - 

z,t  z,l 

i= 1 


(33) 


In  summary  then,  the  K  pairs  of  corresponding  fp’s  f'k  <=>  //“'  are  determined  from  which  the 
linear  2K  set  of  constraints  are  determined  based  on  the  parameters  of  {St,Ct}  .  The  LS 
pseudo-inverse  of  the  constraint  matrix  is  determined  resulting  in  the  LS  estimate  of  {S),C)} 
from  which  the  differential  azimuth  rotation  of  A  a, ,  is  determined.  From  this  the  accumulated 
azimuth  rotation  from  the  first  image  is  determined  as  az  t . 


3-5.  General  affine  transformations 

Next  consider  that  the  camera  is  rotated  and  translated  simultaneously  by  arbitrary  differential 
amounts.  This  is  represented  by  the  affine  transformation  of  the  motion  of  the  kth  feature  point 

from  the  (r-l),/7to  the  t,h  frame  is  described  by  the  affine  transformation.  Instead  of  stating  the 
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affine  transformation  directly  we  start  with  the  perspective  transformation.  Begin  by  defining  the 
perspective  transformation  between  the  static  world  fp’s  as 


.  ,  (34) 

f  =  R,(/  -X) 

Then  we  have 

fW  =  +  T/-!  (35) 
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(37) 


which  is  the  incremental  rotation  from  the  camera  coordinate  in  the  t-1  frame  to  the  world  frame 
and  then  back  to  the  camera  coordinate  in  the  tth  frame.  (!)  -TM)  is  the  incremental 

translation  between  the  t-1  and  t  frames  relative  to  the  world  coordinates.  AsR,_j  rotates  this 
vector  into  the  camera  frame  at  the  (t-1)th  frame,  R;  l  (1)  -TM)  is  the  incremental  translation 
vector  relative  to  the  (t-1)th  frame.  Hence 

i^/  i  (T,  —  T(_j )  -  R,  (Tt  —  T,  j )  (38) 

is  the  translation  vector  relative  to  the  tth  frame.  Hence 
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Now  the  perspective  transformation 
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From  this  the  incremental  affine  transformation  applied  to  the  fp’s  becomes 
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(40) 


From  which  the  differential  affine  transformation  is  defined  as 

C,  S,  Ax, 
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Now  consider  the  cascaded  transformations  of  the  first  two  frames.  We  have 
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Note  that 


C2Cj  -  S2S{  =  cos(AaZj  +  A az2 ) 
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From  this  we  can  define  the  accumulated  affine  transformation  as 
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(44) 
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with  the  components  related  to  the  rotation  and  translation  of  the  camera  relative  to  the  world 
reference  at  the  tth  frame  such  that 
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From  this  we  have  the  accumulated  affine  transformation  given  as  a  recursive  relation  as 


A,=AA,Am  (46) 

There  are  two  approaches  to  estimating  the  differential  affine  transformation  coefficients.  The 
first  is  to  consider  the  2K  constraints  based  on 
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That  is  we  form  for  each  fp  a  pair  of  constraints  as 
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Due  to  changes  in  camera  height,  distortions  and  sources  of  measurement  noise,  in  practice 
C2  +  S2  ^  1 .  This  discrepancy  is  removed  by  the  atan2()  processing  in  regards  to  determining 

the  angle  change.  The  deviation  of  C2  +S2  from  1  provides  an  estimate  of  the  change  in 
height  relative  to  the  current  height. 

As  an  alternative  to  the  LS  solution  of  [Ct,St,Axt,Ayt]  we  can  solve  for  the  six  affine 
coefficients  directly  as  given  by  the  fp  update  of 
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Note  that  there  are  6  unknowns  whereas  previously  there  were  4  unknowns.  However,  the  extra 
degrees  of  freedom  are  not  meaningful  as  there  are  additional  constraints  implied  by  the 
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{Ct,St,Axt,Ayt}  system.  Note  that  the  general  affine  transformation  can  accommodate  rotation 
about  any  arbitrary  center  point  whereas  the  {Cl,St,Axl,Ayl}  system  assumes  rotation  to  be 
fixed  about  the  camera  center.  Also  the  scaling  is  the  same  in  the  xe  and  ye  directions  which 

accounts  for  the  additional  DOF.  To  use  this  system  it  is  necessary  to  account  for  the  rotation 
center  and  the  different  scaling  factors  when  mapping  back  to  the  ground  reference  system. 

The  coefficients  of  {a,b,c,d,e,f}  are  solved  for  with  two  sets  of  equations  as  in  a  previous 

chapter  as 
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3.5  3D  Perspective  ego-motion  based  on  a  2D  marker  and  unknown 
fp's 

Next  consider  that  the  camera  undergoes  6  DOF  motion  change  while  viewing  a  known  2D 
poster  marker  on  a  wall  as  shown  in  Figure  1 1 .  Here  the  marker  is  in  the  Zw  =  0  plane  at  a 

known  location.  Next  consider  the  perspective  mapping  of  a  2D  pattern  on  the  Zw=  0  plane  in 

the  world  coordinates  that  is  imaged  in  the  camera  by  projection  onto  the  image  plane.  Another 
illustration  is  given  in  Figure  3-12  where  a  set  of  fp’s  is  used  to  provide  continuous  3D 
egomotion.  The  concept  is  that  there  will  always  be  several  fp’s  in  the  camera  FOV  which  are 
sequentially  self  calibrating.  That  is  as  a  new  fp  enters  the  FOV  the  existing  fp’s  will  be  used  to 
calibrate  the  position  of  the  new  fp’s.  This  will  give  us  the  perspective  coefficients  which  are 
related  to  {R,T} .  As  Zw  =0we  can  write 
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(52) 


X" 

X" 

RT 

Y* 

=  H 

Yw 

1 

1 

where  H  is  a  3x3  matrix  with  elements  defined  as 


Figure  3-11  Square  poster  viewed  from  a  camera  in  the  image  plane  can  be  restored  to  the 
undistorted  square  shape  based  on  the  perspective  transformation. 
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R  -  rotation  between  world  and  camera  reference 


T  -  translation  vector  between  world  and  camera  reference 


trajectory  of  camera 


Figure  3-12  Concept  of  3D  ego-motion  based  on  viewing  the  floor  surface  with  various  fp ’s  in 
the  current  camera  FOV 


Note  that  the  first  two  columns  of  H  are  the  first  two  columns  of  R  and  the  third  column  of  H  is 
-RT  as  before.  Every  fp  of  the  marker  that  is  distinguishable  provides  two  constraints  from 
which  the  9  components  of  Hcan  be  determined.  We  have 

X  =  x  ,z  HuXw  +  HuYw+Hl3 
cl  c  H31XW+H32YW+H33 
Y/z  H^+H^+H,  3 

c/  c  H3lXw  +  H32Yw  +  H33 

which  is  rearranged  as 

-HnXw-HuYw-Ha  +H3lxXw  +  H32xYw+H33x  =  0 
-H2lXw  -  H22Yw  -  H23  +  HnyXw  +  H32yYw  +  H33y  =  0 

This  results  in  a  pair  of  constraints  expressed  as 
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(56) 


uxb  =  0 
Uyb  =  0 

where 


b  = 

\Hn  Hn 

h21  H22 

^23 

h31 

H,2 

Ux 

II 

1 

* 

1 

0  0  0 

*YW 

A 

Uy 

II 

o 

o 

o 

1 

-K  -i 

yYw 

y\ 

T 


For  a  set  of  K  fp’s  that  can  be  identified  from  the  marker  we  have 


(57) 


u 


X,1 


u 


y.i 


u 


\,K 


u 


y  ,k 


b  = 


(58) 


We  then  apply  the  singular  value  decomposition  (svd)  solution  to  determine  the  components  of 
b  based  on  using  the  right  singular  vector  corresponding  to  the  smallest  singular  value.  Once 
the  components  of  H  is  determined  then  we  can  map  these  into  the  set  of  variables 
{ax,ay,az,xT,yT,zT,s}ftom  which  the  {R,T}  matrices  can  be  determined.  The  additional 

variable  of  s  is  a  scaling  factor  that  is  necessary  as  Hwill  generally  have  an  arbitrary  scaling 
associated  with  it.  Hence  s  is  a  nuisance  parameter  but  necessary  in  order  to  make  the 

mapping  of  H— >•{ ax,ay,az,xT,yT,zT}  possible.  The  following  Matlab  routine  shows  the 
procedure  of  extracting  {ax,ay,az,xT,yT,zT}  from  H. 

First  we  generate  an  H  matrix  to  work  with. 

%  Hmat 

function  H  =  Hmat (ax,ay,az,T,s) 


%  generate  the  RT  matrix 
Cx  =  cos (ax) ;Sx  =  sin  (ax) ; 
Cy  =  cos(ay);Sy  =  sin (ay) ; 
Cz  =  cos (az) ;Sz  =  sin(az); 


Rx  =  [ [ 1 , 0 , 0 ] ;  [ 0 , Cx,  Sx] ;  [ 0 , -Sx,  Cx] ]  ; 
Ry  =  [ [Cy  0  Sy] ; [0  1  0] ; [-Sy  0  Cy] ] ; 
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Rz  =  [ [Cz  Sz  0 ] ;  [  -Sx  Cz  0 ]  ;  [ 0  0  1]  ]  ; 

R  =  Rz*Ry*Rx; 

H  =  s* [R ( :  ,  1 : 2)  ,  -R*T]  ; 

The  function  to  determine  the  error  vector  used  in  the  nonlinear  equation  solver  is  given  as 

function  F=H2RTfun (H, x) 
ax  =  x ( 1 ) ; 
ay  =  x (2) ; 
a  z  =  x  ( 3 )  ; 

T  =  x  (4  :  6)  ; 
s  =  x  ( 7  )  ; 

e  =  H-Hmat(ax,ay,az,T,s)  ; 

F=zeros (9,1) ; 

F ( : )  =  e; 

Finally  the  call  to  fsolve(),  Matlab’s  generic  equation  solver  is 

options=optimset (' Display iter ')  ;  %  Option  to  display  output 

[aTc, fval , exit flag, output, j  acobian]  =  f solve ( @ (x) H2RTfun (H, x) , aT, options ) ; 

The  fsolve()  routine  has  no  issue  with  quickly  converging  on  the  correct  values  of  the 
parameters. 

Next  consider  the  more  general  problem  of  the  markers  unknown  to  the  camera.  However,  to 
simplify  the  scenario  it  will  be  assumed  that  the  perspective  transformation  is  known  at  the  initial 
point  where  the  camera  begins  it’s  trajectory.  At  this  point  it  then  begins  to  observe  unknown 
markers.  It  is  also  assumed  that  the  initial  known  template  and  fp’s  are  all  on  the  Zw  =  0  plane. 

Consider  that  the  camera  undergoes  an  arbitrary  translation  and  rotation  relative  to  the  world 
reference  such  that  at  the  time  that  the  tth  frame  is  available  from  the  camera  the  perspective 

rotation  matrix  and  translation  vector  are  denoted  as  {R,,T/} .  The  objective  is  that  given 

{Rm,T(_,}  for  the  (t-1  )th  frame  then  determine  {R,,Tj  for  the  tth  frame  based  on  the  movement 

of  the  fp’s  in  the  image  plane  as  recorded  by  the  camera.  As  derived  before,  the  transformation 
between  the  two  camera  positions  is  given  as 


(59) 
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where  Pj  1  and  Pc'  are  the  position  vectors  of  the  fp  in  the  camera  reference  in  the  (t-1)th  and  tth 
frame.  Define 


(60) 

as  the  rotation  matrix  between  the  two  camera  coordinates  from  frame  t-1  to  frame  t.  That  is 


x'nx*-1 

X‘QVf' 

x‘cz; 

7S 

>  ^ 

II 

Yc'rx‘ 1 

Yc‘[]YcM 

Yc‘nz; 

z*ax‘-‘ 

z^;-1 

z*m: 

Also  define 

11=Rm(T,- 

Tm) 

(61) 


(62) 


as  the  translation  vector  between  the  two  camera  positions  relative  to  the  camera  frame  at  t-1. 
Then  the  transformation  is  given  as 


(63) 


Note  that  Pc'  1  is  known  as  {Rm,T(_J  as  well  as  the  image  coordinates  of  the  fp  from  the  (t-1)th 
frame  are  known.  We  can  demonstrate  this  as  we  have  the  following  equations 


c,x 

( 

p 

1  W,x 

p,-\ 

c,y 

p,-\ 

=  R1 

p 

W,y 

0 

-T 

A/-i 

c,z 

\ 

J 

v-1 =<;7C1 
y'-'  =  /  p:: 


(64) 


which  represents  5  equations  and  5  unknowns  as  [P^ ^ y  ' ,PW X,PW y}  •  s  is  a  known 
scaling  factor. 

Define  H  as  a  3x4  matrix  as 
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h=[r;  -r;t']= 


Kx  Ki  Ki  K 


h2 1  h 


2,2 

h3l  h32  h3  3 


h2  3  h 


2,4 

^3,4 


(65) 


The  non-homogeneous  image  coordinates  of  the  tth  frame  are 

P' 


x  =  s  c,x 


P' 

t  c,y 

yc=s- 


(66) 


P 


Let  hi=\[hil  hj2  hi3  /*.4]  denote  the  ith  row  of  H  and  let  /  =  ^Pc'x1  pj^1  P'J  lj  denote 
the  homogeneous  fp  coordinates  such  that 


Kfc  =  Kf 
KM  =  Kf 

This  results  in  a  pair  of  homogeneous  equations  for  each  fp  expressed  as 


/  0  0  0  0  —Jx‘c 

0  0  0  0  /  -Jy‘c 


(67) 


iK" 

"0" 

hT, 

— 

z 

0 

\hl_ 

(68) 


From  this  a  pair  of  constraints  are  formed  for  each  fp  correspondence  which  can  be  expressed 
as  a  homogeneous  solution  of  the  unknown  values  of  htj .  First  a  1x12vector  of  htj  denoted  as 

h  is  defined  as 


*  =  [*,  K  h,  f 

=  [K  Ki  Ki  Ka  K\  Ki  K3  Ka  K\  Ki  Ki  Ka\ 

Set  up  the  matrix  component  for  the  kth  fp  as 


(69) 


- 1 

o 

o 

o 

o 

1 

_ 1 

1 

_ 1 

1  o 

o 

o 

o 

1 

>> 

1 _ 

1 

A 

1 _ 

(70) 
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where/;  =[P;-\ 


Pet  l] 


Note  we  start  with  the  assumption  that  we  have  P .  Then  from  the  corresponding  fp’s  that  we 

observe  in  frame  t  we  can  solve  for  the  matrix  H .  From  the  matrix  H  we  get  and  .  Then 
we  can  solve  for  the  rotation  matrix  as 


R/=R'aR,-i  (71) 

and 

T,=RlX+T,_t  (72) 

Then  we  have 

P,'=Ri(pr'-T4)  (73) 

and  we  are  set  for  another  iteration. 

The  overall  algorithm  can  be  described  as  follows.  First  the  known  marker  is  observed  in  the 
camera  at  the  known  location  for  which  we  fix  the  world  reference.  From  this  the  initial  {R0,T0} 
is  established.  Next  the  camera  also  observes  some  fp’s  of  opportunity  at  unknown  location 
except  for  the  assumption  that  the  fp’s  are  located  also  on  the  Zw  =  0  plane.  Through  multiple 
consecutive  frames,  while  both  the  known  marker  and  fp’s  are  in  the  camera  FOV,  the  location 
of  the  fp’s  is  determined  in  the  world  reference  in  the  Zw  =  0  plane.  As  the  camera  is  moved, 

the  known  marker  will  eventually  drift  out  of  the  FOV  and  we  have  only  the  fp’s  of  opportunity. 
However,  as  the  location  of  these  have  been  established,  then  we  can  determine  the 
perspective  transformation  to  the  camera  and  also  {R,,TJ  .  As  the  trajectory  progresses,  new 
fp’s  are  going  to  enter  the  FOV  with  positions  estimated  based  on  the  older  set  of  fp’s  and 
hence  the  process  of  extracting  {R(,T(}  can  continue  indefinitely  as  long  as  the  fp’s  are  all  on 

the  Zw  =0  plane. 

In  a  practical  scenario  there  will  be  limitations  and  complications  with  this  method.  Due  to  noise 
and  lens  distortion  of  the  camera  there  will  be  a  drift  issue.  That  is  as  the  trajectory  progresses 
and  new  fp’s  are  added  to  the  mix  then  the  errors  will  accumulate.  Also  the  camera  will  pick  up 
fp’s  that  are  not  on  theZw  =  0  plane.  These  will  have  to  be  processed  based  on  determining  the 
plane  in  the  world  reference  where  they  reside. 
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3.7  Alternative  3D  Method 


While  the  previous  method  works  in  general  3D  as  will  be  shown  experimentally,  it  is  subject  to 
drift  as  it  is  necessary  to  base  {R, ,  T, }  on  the  previous  step  of  {R,_, ,  T,_, } .  Additionally  it  is 


necessary  that  the  initial  pose  of  the  camera  is  known  as  eventually  R0  is  assumed  to  be 

known.  Generally  3D  egomotion  of  this  type  is  based  on  epipolar  geometry.  However,  to 
determine  the  essential  matrix  for  the  epipolar  update,  it  is  necessary  to  have  fp’s  that  are  not 
on  a  plane.  At  the  same  time  the  powerful  constraint  of  the  fp’s  residing  on  the  plane  is  not 
utilized.  Hence  another  approach  is  considered  which  specifically  takes  advantage  of  the 
Zw=  0  constraint.  Consider  that  we  have  two  images  from  two  camera  positions.  For  a  given 


fp  we  have  the  position  vector  from  camera  1  given  as 

7 


Xx=R, 


fy 

0 


+z 


where  is  the  rotation  from  the  planar  world  to  the  camera  in  position  1  and7j  is  the  translation 

from  the  world  to  camera  1  referenced  to  the  camera  coordinate.  We  have  the  image  plane 
coordinates  of  the  camera  in  position  1  such  that  the  observable  is 


x,  = 


i  -  y  ~XU 


(28) 


where  s  is  the  scaling  factor.  Likewise  we  have  for  camera  position  2 


X2  =R2 


7, 

fy 


+i ; 


where  R2  is  the  rotation  from  the  planar  world  to  the  camera  in  position  2  and  T2  is  the 

translation  from  the  world  to  camera  2  referenced  to  the  camera  2  coordinates.  Image 
coordinates  are  given  as 
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(29) 


y  i  — 


x. 


2,z 


Next  consider  the  transformation  of  camera  in  position  2  relative  to  camera  position  1  given  as 


X2=RrXl+Tr 

7" 

=  RA  fy  +  KT\  +  Tr 


(30) 


0 

is  the  rotation  from  camera  frame  1  to  camera  frame  2  in  frame  2.  Tr  is  the  translation  from 
camera  1  to  camera  2  in  the  frame  of  camera  2.  That  is  relation 


(31) 


where  Xk .  represents  the  ith  basis  vector  of  camera  position  k.  Also 

Tr=ox-o2 

where  is  the  camera  center  of  the  kth  position  relative  to  the  camera  position  2  frame.  Now 
we  define  d  is  the  perpendicular  distance  from  the  plane  to  ox  (camera  position  1 )  and  N  is  the 

normal  to  the  plane  of  Zw  =  0  relative  to  camera  frame  1  such  that  d  =  NT Xx .  Note  that  there  is 
an  ambiguity  as  the  normal  to  the  world  plane  N  can  have  two  ambiguous  directions.  This  is 
resolved  by  stipulating  that  d  >  0 .  A  trick  is  to  set  1  =  NTXx/d  such  that 


x2=Rrxl+i-Tr 
„  „  (  NTX , 


(32) 
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(Note  the  order  of  the  vectors  can  be  changed  as  they  are  dotted  to  a  scalar.)  Then  we  have  a 
matrix  that  relates  these  two  positions  as 


X2  =  HXx 

1  r  (33) 

HURr+-TrNT 
d 

As  observed,  we  have  the  ratio  T/d which  is  dimensionless.  Hence  we  cannot  determine  both 

the  height  of  the  camera  and  the  translation  simultaneously  but  only  the  ratio.  Next  scale  the 
position  vectors  by  normalizing  to  form  the  2D  homogeneous  vectors  as 


=  Alxl  = 


x2  —  ^2^2  — 


7i 

1 


72 

1 


We  can  then  write  X2  =HX \  relating  the  2D  homogeneous  vectors  as 


(34) 


x2  =  AHx  j  (35) 

Note  that  the  scaling  constant  A  is  irrelevant  as  it  cannot  be  determined  anyways  since  H  will 
only  be  known  to  within  a  scaling  constant.  Now  consider  that  x2  xx2  =  x2x2  =  0  such  that 


x2Hx  j  =  0 


(36) 


which  is  a  set  of  homogeneous  relations.  Note  that  x2  can  be  expressed  as  a  skew  symmetric 
matrix  operator  as 

0  -z  y 


x  = 


z 


0  -x 


-y  x  o  _ 

To  solve  for  the  elements  of  H  within  a  scaling  constant  the  following  procedure  is  used 
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(37) 


aJ  =  x{  0 

A  =  [a1  -  aKJ 

Solve  the  SVD  of 


AHS  =  0 


(38) 


where  Hs  is  the  stacked  version  of  H  corresponding  to  the  right  singular  vector  with  the  smallest 
singular  value.  The  second  singular  value  of  H  should  be  1  which  gives  a  way  of  normalizing  H. 

1  T  1 

After  determining  //then  H  =  Rr+  —  TrN  determine  i?rand  —  Tr  .  Note  that  we  can  only  get  a 

d  d 

scaled  version  of  Tjd  such  that  the  height  is  required  or  some  measure  of  distance  for 
calibration.  Another  complication  is  that  N  is  in  the  camera  position  1  frame. 

Example  -  Consider  the  simplest  example  of  a  simple  translation  of  tx  to  the  right  where  the 
plane  of  the  fp’s  is  a  distance  d  from  the  camera  positions  and  the  normal  to  the  plane  is 
iV  =  [0  0  if .  (Note  that  N  *  [0  0  -if  as  NTXl  =d).  We  have  Tr  =  \-tx  0  of  and 


Rr=I3.  Hence, 


H  =  Rr+-TN7 
a 


'1 

- 

tjd 

= 

1 

+ 

0 

1 

0 

1 

o 

1 

_ i 

= 

0  1 

0 

0  0 

1 

i] 


From  x2  =/l//x1we  have 


x2 

'1 

0 

~-tJd 

*1 

1 

>< 

1 

vT* 

1$: 

72 

=  A 

0 

1 

0 

7i 

=  A 

7i 

1 

0 

0 

1 

1 

1 

which  makes  sense. 

Next  consider  a  simulation  of  determining  H  arbitrary  rotation  angles  and  translation.  Again  we 
consider  a  set  of  random  feature  points  distributed  in  the  world  plane  around  the  origin.  Note 
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that  care  has  to  be  taken  that  the  translation  vector  is  referenced  to  the  correct  coordinate 
system. 

First  define  the  parameters,  d  is  the  height  of  the  camera  1  position,  (tx,ty,tz)  are  the 
components  of  the  translation  vector  in  the  coordinates  of  position  1,  sc  is  the  scaling  of  the 
camera  to  pixel  index,  and  (ax,ay,az)  are  the  rotation  angles  of  position  2  relative  to  position  1 . 
Provided  that  sc  is  the  same  in  the  x  and  y  directions  it  has  not  impact  on  the  estimation  of  H. 

d  =  10/  %  height  of  camera  for  position  1 

tx  =  1 ; ty=0 . 5; tz=0 ;  %  relative  motion  between  positions  1  and  2 

ax  =  0.1; ay  =  0.2;az  =  0.1; 
sc  =  1;  %  scaling  of  camera 
K  =  10;  %  number  of  random  points 

Next  we  generate  the  relative  rotation  and  translation  matrices  denoted  as  Rr  and  Tr 
respectively. 

Cx  =  cos (ax)  ;Sx  =  sin  (ax); 

Cy  =  cos(ay);Sy  =  sin  (ay); 

Cz  =  cos(az) ;Sz  =  sin(az); 

Rx  =  [ [ 1 ,  0 ,  0 ] ; [ 0 ,  Cx, Sx] ; [ 0 , -Sx, Cx] ] ; 

Ry  =  [ [Cy  0  Sy] ; [0  1  0] ; [-Sy  0  Cy] ] ; 

Rz  =  [[Cz  Sz  0];[-Sz  Cz  0] ; [0  0  1]  ]  ; 

Rr  =  Rz*Ry*Rx; 

Tr  =  [ tx; ty; tz ] ; 

The  position  1  is  given  the  translation  of  T_1  and  R_1  relative  to  the  position  1  frame 

T_1=R_1  *  [ 0 ; 0 ; d ] ; 

R_1  =eye (3 ) ; 

H_1  =  [[sc  0  0] ; [0  sc  0] ; [0  0  1]]  *  [R_l ( : , 1 : 2 ) , R_1*T_1 ] ; 

The  position  2  is  given  the  translation  of  T_2  and  R_2  relative  to  the  position  2  frame 

T_2  =  R_2* (T_l  +  Tr)  ; 

R_2  =Rr*R_l; 

H_2  =  [[sc  0  0 ] ; [ 0  sc  0 ] ; [ 0  0  1]]  *  [R_2 ( : , 1 : 2 ) , T_2 ] ; 

Generate  a  set  of  random  feature  points  and  map  these  into  position  1  and  position  2  images. 

px  =  3*randn (K, 1 ) ; 
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py  =  3*randn (K, 1 ) ; 


%  Determine  the  images  as  seen  in  camera  1  and  camera  2  positions 
xl  =  zeros (K, 1 ); x2  =  zeros (K, 1 ); yl  =  zeros (K, 1 ); y2  =  zeros (K, 1) ; 
for  k=l:K 

gl  =  H_l* [px (k) ;py (k)  ;  1]  ; 
xl (k)  =  sc*gl (1) /gl (3) ; 
yl (k)  =  sc*gl (2) /gl  (3) ; 
gl  =  H_2* [px (k) ;py (k) ; 1] ; 
x2(k)  =  sc*gl (1) /gl (3) ; 
y2 (k)  =  sc*gl (2) /gl  (3) ; 

end 


Generate  the  a1  =  x(  ®x]2  components  and  put  them  into  the  matrix  A  =  ^al  •••  aK~^ 

A  =  zeros (3*K, 9) ; 
for  k=l:K 

a=kron ( [xl ( k) ; yl  ( k) ; 1 ] ,  [  [0  -1  y2(k)];[l  0  -x2 (k) ] ;  [-y2 (k)  x2 (k)  0]]) 

A ( ( ( k-1 ) *3+1) : k^3 , : ) =a ’ / 

end 


Solve  the  SVD  of 


AHS  =  0 


as 

[V,  L,  U]  =  svd  (A)  ; 

H  =  zeros (3) ; 

H  ( :  )  =  U(:,9)  ; 
gl  =  svd (H) ; 

H  =  sign (H  (3 , 3)  ) *H/gl (2 ) 


Extraction  ofi?rand  —Tr 

d 


Use  fsolve()  to  determine  i?rand 


d 


Tr  from  H  given  [d,  N) . 


options=optimset (' Display ' ,  ’ iter ’) ;  %  Option  to  display  output 

%  initial  guess  set  to  first  location 
x  =  [0; 0; 0; Tr; 1]  ; 

[x, fval , exit flag, output, j  acobian]  = 

f solve ( @ (xx) Hr2RrTrf un (H, xx,  d,  [  0 ;  0 ;  1  ]  )  ,  x, options )  ; 

D=x; 

a_e  s  t  =  x ( 1 : 3 ) 

Tr_est  =  R_2'*x(4:6) 

function  F=Hr2RrTrf un (Hr , x, d, N) 
ax  =  x(l);ay  =  x  (2)  ;az  =  x(3); 

Tr  =  x(4:6);s  =  x  ( 7 )  ; 
e  =  Hr-Hrmat (ax, ay, az , Tr , s , d, N) ; 

F=zeros (9,1) ; 

F ( : )  =  e; 

function  H  =  Hrmat (ax, ay, az , Tr , s , d, N) 

%  d  is  the  perpendicular  distance  from  ZW=0  plane  to  position  1 
%  N  is  the  normal  vector  of  the  plane  ZW=0  in  pos  1  frame 
%  generate  the  RT  matrix 
Cx  =  cos (ax) ;Sx  =  sin (ax); 

Cy  =  cos(ay);Sy  =  sin (ay); 

Cz  =  cos (az ) ; Sz  =  sin(az); 

Rx  =  [  [1, 0, 0 ] ;  [ 0, Cx, Sx] ;  [ 0, -Sx, Cx] ] ; 

Ry  =  [ [Cy  0  Sy] ;  [0  1  0] ;  [-Sy  0  Cy]  ]  ; 

Rz  =  [ [Cz  Sz  0] ;  [-Sx  Cz  0] ;  [0  0  1] ]  ; 

Rr  =  Rz*Ry*Rx; 

H  =  s* (Rr+Tr*N ’ / d) ; 

3.8  Experimental  results  of  3D  ego-motion  based  on  Tracking  fiducial  markers 

In  this  section  some  of  the  experimental  results  will  be  given  that  are  based  on  the  tracking  of 
fiducial  markers.  The  concept  is  illustrated  in  a  smaller  scale  in  the  photo.  The  handheld 
naviagation  device  (HND)  observes  some  fiducial  markers  in  its  FOV.  It  can  generate  the 
perspective  matrix  of  the  relative  3D  position  of  the  camera  with  respect  to  the  fiducial  marker. 
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Figure  3.12  Figure  depicting  the  camera  and  world  coordinate  system 


As  the  camera  is  moved,  the  initial  fiducial  will  leave  the  FOV  but  other  fp’s  on  the  same  surface 
will  be  observed  as  illustrated  in  Figure  3.13.  The  perspective  transformation  determined  based 
on  the  first  fiducial  can  be  used  to  locate  the  fp’s  of  opportunity  seen  in  the  FOV  at  the  same 
time.  The  fp’s  located  are  then  used  as  an  equivalent  fiducial  marker  and  so  forth.  Eventually 
another  calibrated  fiducial  marker  known  to  the  camera  will  emerge  in  the  FOV  and  the  drift  can 
be  corrected. 
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Figure  3.13  Camera  is  translating  illustrating  a  continual  selection  of  fp ’s  in  the  FOV 

Figure  14  shows  an  image  of  a  fiducial  marker  that  is  seen  in  perspective  and  corrected  after 
the  transformation  has  been  determined  and  applied.  The  projective  transformation  that  is 
determined  is  then  used  to  infer  the  ego-motion  of  the  camera  itself  relative  to  the  fiducial 
marker. 
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Figure  3.14  Top  image  is  fiducial  in  perspective  from  which  the  projective  transformation  is 
determined  and  applied  resulting  in  the  transformed  orthographic  view 

The  setup  for  the  trial  with  the  sequence  of  fiducial  markers  and  posed  fp’s  of  opportunity  is 
shown  in  the  photo  of  Figure  3.15. 


Figure  3.15  Experimental  setup  for  the  sequence  of fiducial  markers  and fp  ’s  of  opportunity 

The  results  of  this  trial  is  shown  in  Figure  3.16  for  one  of  the  runs.  The  results  show  some  drift 
that  occurs  as  well  as  the  issue  with  the  small  discontinuities  in  the  solution  as  the  fiducial 
markers  enter  and  leave  the  FOV.  The  drift  is  an  annoyance  but  it  is  not  large  and  hence  is 
manageable.  When  other  markers  are  added,  the  drift  will  correct  itself.  Also  the  C V  method  is 
not  intended  as  the  only  source  of  navigation  measurements.  Wireless  sources,  as  explained  in 
the  previous  chapters  are  less  accurate  over  the  short  trajectory  but  do  not  suffer  from  the  same 
drift  issue. 
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Transition  in  FOV  from  FM  lo  FMO  -  Top  view 


Outcome  of  a  longer  length  trial  is  shown  in  Figure  3.17  in  the  same  hallway  setting.  Drift  is  still 
present  but  the  positioning  from  the  3D  projective  is  useably  accurate. 
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Figure  3.17  Experimental  outcome 


A  least  squares  fit  is  applied  to  the  trial  results  of  Fig. 3-1 7  and  are  given  in  Figure  3-18.  The 
least  squares  is  fit  to  each  segment  with  the  same  set  of  fiducial  markers  and  features.  The 
overall  least  squares  shows  very  little  drift.  As  the  camera  was  mounted  on  a  cart  of  constant 
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height  from  the  floor  surface,  the  height  should  ideally  be  constant  which  it  is  approximately  with 
very  little  drift  with  respect  to  this  variable. 
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Figure  3.1 7  Experimental  outcome  after  least  squares  line  fit 


Having  the  camera  determine  the  various  projective  transformations  the  observations  after  each 
frame  will  be  entered  into  the  EKF  based  on  a  linear  model  for  the  camera  trajectory  as  outlined 
in  the  previous  chapter.  The  processing  block  diagram  is  given  in  Figure  3.18. 
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Figure  3-18  Kalman  filter  processing  used  to  determine  the  least  squares  trajectory  of  the  HND 

3.9  Extracting  Lines  with  Hough  Transforms 

In  addition  to  the  fp’s  of  the  wall  or  floor  surface  and  fiducial  markers  it  is  convenient  to  use 
rectangular  structures  such  as  windows,  doors  and  picture  frames  as  well  as  lines  that  exist  in 
every  room.  As  projective  transformations  of  lines  in  3D  are  lines  in  the  image  plane  such 
structures  are  preserved.  Figure  3-19  illustrates  some  of  the  features  and  lines  that  are 
available  in  a  3D  scenario. 
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Figure  3-19  Typical  indoor  setting  has  extensive  instances  of  lines  and feature  points. 


Geometric  structures  such  as  lines,  rectangles  and  circles  can  be  found  in  a  camera  image 
using  the  general  Hough  transform  It  is  essentially  a  transform  that  determines  the  likelihood  of 
the  geometric  construct  appearing  in  the  image  over  the  range  of  the  unknown  variables.  For 
instance,  consider  a  line  which  is  defined  by  two  parameters,  the  slope  angle  and  skew  distance 
from  the  origin.  The  Hough  transform  would  map  the  image  in  x-y  to  the  two  dimensional 

uknown  parameter  space  of  slope  angle#  and  skew  distance p .  A  transform  point  of  (p,9) 

then  is  an  indication  of  the  likelihood  that  a  line  with  the  parameters  of  (p,9)  exists  in  the 

image.  For  a  Hough  ellipse  transform,  there  are  five  parameters.  Usually  the  transform  is  a 
construction  of  two  of  these  parameters  such  that  the  other  three  are  given.  For  instance  the 
two  parameters  could  be  the  center  of  the  ellipse  with  the  major  and  minor  axis  given  as  well  as 
the  orientation  angle.  This  is  better  shown  with  worked  examples  of  the  line  and  the  circle. 

As  shown  in  the  figure  below,  any  line  segment  that  is  extended  into  an  infinite  line  can  be 
described  by  two  variables: 
p  -  skew  distance  from  the  origin  to  the  line 

9  -  angle  of  the  skew  line 

For  each  point  in  the  intensity  field  plot  a  contour  of  all  the  potential  { p,0 }  combinations  or  lines 

that  the  intensity  point  could  belong  to.  This  line  is  weighted  by  the  intensity  of  the  point.  In  this 
way  when  sweeping  over  an  intensity  field,  the  points  corresponding  to  a  bright  line  of 

parameters  { p,0 }  will  add  up  to  form  a  peak  in  the  2D  intensity  field  of  [p,9].  This  intensity 

plot  is  then  thresholded  to  yield  the  estimate  of  the  { p,0 }  of  the  line  segment.  Note  that  the 
longer  and  more  intense  this  line  is,  the  higher  the  peak  will  be. 
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Figure  6-10.  A  point  (je0J  y0)  in  the  image  plane  (panel  a)  implies  many  lines  each  parameterized  by 
a  different  p  and  8  (panel  b);  these  lines  each  imply  points  in  the  ( p ,  8)  plane,  which  taken  together 
form  a  curve  of  characteristic  shape  (panel  c) 


Figure  3-20  Hough  line  detection  (taken  from  Learning  OpenCV  Bradlie) 


An  example  of  using  the  Hough  transform  for  lines  based  on  the  image  in  Figure  3-21  is  given. 
First  the  image  is  converted  to  grey  scale  and  then  the  Canny  edge  detector  is  applied  (taken 
from  Matlab,  Image  processing  toolbox  user  guide. 


Figure  3-21  Image  of  gantry  section 
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Figure  3-22  Edge  detection  of  gantry  section 

Next  take  the  Hough  line  transform  of  the  edge  image  which  is  shown  in  Figure  15. 

[H, T, R]  =  hough (BW) ; 

figure ( 4 ) ; imshow (H, [ ] ,  ' XData ' ,  T ,  ' YData ' , R, . . . 

' InitialMagnif ication ’ ,  ' fit ' ) ; 
xlabel ( ’ \theta ’ ) ,  ylabel ( ' \rho ’ ) ; 
axis  on,  axis  normal; 


H  is  the  likelihood  as  a  function  of  p  and  6 .  It  can  be  plotted  as  a  mesh  or  contour  which 
sometimes  provides  greater  clarity.  The  highest  peaks  of  H  which  correspond  to  the  most 
probable  lines  are  given  by  houghpeaks  which  can  be  added  to  the  image  using 

hold  on; 

P  =  houghpeaks (H, 5 ,  ' threshold ceil ( 0 . 3*max (H (:)))) ; 

x  =  T(P(:,2)  )  ;  y  =  R(P(:,1)); 

plot (x,  y,  's’,  'color ' ,  'white ' )  ; hold  of f ; 
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Figure  3-23  Hough  transform  of  edge  image  with  highlighted  peaks  around  0  =  50' 


Now  find  the  lines  and  map  them  onto  the  figure 

%  Find  lines  and  plot  them 

lines  =  houghlines ( Ic, T, R, P,  ’ FillGap ’ ,  5 ,  ’ MinLength ’ , 7 ) ; 

figure,  imshow(I),  hold  on 

max_len  =  0/ 

for  k  =  1 : length ( lines ) 

xy  =  [ lines ( k)  . pointl ;  lines ( k)  . point2 ] ; 

plot (xy ( :  ,  1 )  ,  xy ( :  ,  2  )  ,  ’ LineWidth ’  ,  2 ,  ’Color ' ,  ’ green ’ )  ; 

%  Plot  beginnings  and  ends  of  lines 

plot (xy (1,1) , xy (1,2) ,  ’ x ’ ,  ’ LineWidth ’ , 2 ,  ’Color ’  ,  ’yellow ' )  ; 
plot (xy (2,1) , xy (2,2)  ,  ’ x ’ ,  ’ LineWidth ’ , 2 ,  ’Color ’  ,  ’ red ’ )  ; 

%  Determine  the  endpoints  of  the  longest  line  segment 
len  =  norm ( lines ( k) . pointl  -  lines ( k) . point2 ) ; 
if  (  len  >  max_len) 
max_len  =  len; 
xy_long  =  xy; 

end 

end 

%  highlight  the  longest  line  segment 

plot (xy_long ( : , 1 ) , xy_long ( : , 2 ) , ’ LineWidth ’ , 2 , ’ Color ’ , ’ blue ’ ) ; 


Figure  3-24  Hough  lines  mapped  back  into  image  with  starting  and  stopping  positions 

Note  that  there  are  five  lines  but  that  each  line  consists  of  multiple  segments.  Hence  hough  lines 
finds  the  line  segments  within  the  line.  The  segmented  line  is  an  additional  step  beyond  the 
Hough  transform.  Line  segments  are  based  on  the  intensity  of  the  edge  detection  of  the  image 
itself. 
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Figure  3-25  Hough  Lines  applied  to  the  building  image 

An  interesting  development  for  indoor  navigation  is  based  on  tracking  floor  tiles  using  Hough 
transform  and  projective  transformations.  Initially  try  the  features  of  opportunity  in  the  video 
which  is  not  successful  as  seen  in  Figure  3-26. 


Figure  3-26  Tiled floor  surface  with  corners()  applied  to  determine  features  of  opportunity 

There  is  an  obvious  difficulty  in  establishing  the  movement  of  features.  After  some  careful 
manual  adjustment  of  the  threshold  gives  us  some  hope  of  using  a  better  set  of  feature  points 
that  can  be  used  for  differential  movement  as  shown  in  Figure  3-27  where  the  green  and  red 
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fp’s  indicated  are  for  two  consecutive  image  frames.  Clearly  it  is  not  a  robust  solution, 
especially  as  manual  search  for  the  optimum  threshold  is  required. 


Figure  3-27  Tiled  floor  with  fp ’s  of  two  consecutive  frames 

With  Gaussian  spatial  filtering  and  tweaking  of  quality  threshold  we  can  get  fp’s  from  the  corners 
of  the  grout  lines.  A  problem  is  that  the  fp’s  jump  around  as  the  metric  for  the  fp’s  has  a  saddle 
point  in  each  grout  line  intersection  making  the  location  of  the  fp  very  sensitive  to  any  noise  or 
neighbourhood  clutter. 


ka 


Figure  3-28  Tiled floor  feature  points  after  Gaussian  filtering  applied 

Next  we  try  the  Hough  line  transform  and  we  get  the  plot  of  the  transform  in  Figure  3-29  and  an 
overlay  spatial  image  as  in  Figure  3-30. 
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Figure  3-29  Hough  line  transform  output  for  tiled floor  image  showing  sets  of peaks  indicating 
the  vertical  and  horizontal  lines 


Figure  3-30  Hough  line  transform  output  for  tiled floor  image 

Next  we  want  to  extract  the  infinite  lines.  This  works  well  if  we  can  assume  that  the  camera  is 
well  calibrated  and  that  the  tiling  is  accurate.  Note  that  the  Hough  line  output  is  P  which  are 
indices  of  T  and  R.  and  T  are 


figure (4),  imshow ( tile_vid ( : , : , 1 ) ) ,  hold  on 
NL  =  length (P); 

Nx  =  640; Ny  =  480; 
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for  j=l:NL 

theta  =  T (P ( j , 2 ) ) *pi/180 ;  roe  =  R(P(j,l)); 

xo  =  roe*cos (theta) ; yo  =  roe*sin (theta) ; xl  =  -sin (theta) ; yl 
pp  =  zeros (4,3) ; 

%  Solve  for  x=l  intercept 
t  =  ( 1-xo) /xl ; 
ya  =  yo+t*yl; 
if  ya>0  &  ya  <  Ny 

PP  ( 1 ,  : )  =  [ 1 , ya , 1 ] ; 

else 

PP(1,  :)  =  [l,ya,0]  ; 

end 

%  Solve  for  x=Nx  intercept 
t  =  (Nx-xo) / xl ; 
yb  =  yo+t*yl; 
if  yb>0  &  yb  <  Ny 

pp  (2,  :  )  =  [Nx,  yb, 1]  ; 

else 

pp  (2 ,  :  )  =  [Nx,  yb, 0 ]  ; 

end 

%  Solve  for  y=l  intercept 
t  =  ( 1-yo) /yl; 
xa  =  xo+t*xl; 
if  xa>0  &  xa  <  Nx 

pp (3, : )  =  [xa, 1,1]; 

else 

PP  (3,  :)  =  [xa,  1 , 0  ]  ; 

end 

%  Solve  for  y=Ny  intercept 
t  =  (Ny-yo) / yl ; 
xb  =  xo+t*xl; 
if  xb>0  &  xb  <  Nx 

pp  (4,  : )  =  [xb, Ny, 1] ; 

else 

PP  ( 4 ,  :)  =  [xb, Ny, 0 ] ; 

end 


%  Plot  the  line 

PPP  =  [ ] ; 

for  j  j  =  1:4 

if  pp  ( j  j , 3 )  ==  1; 

ppp  =  [ PPP ; PP ( j  j , 1 : 2 ) ] ; 

end 

end 

plot (ppp ( : , 1 ) , ppp ( : , 2 ) , ' LineWidth ' , 2 , ' Color ' , ' green ' ) ; 
end 

hold  off; 


cos (theta) ; 


106/143 


Figure  3-31  Hough  line  transform  output  for  tiled floor  image  showing  grid  of  infinite  lines 

With  the  Hough  lines  the  intersection  points  can  be  calculated  for  the  grid  of  tiles.  Note  that 
each  Hough  line  is  an  average  of  all  of  the  points  that  have  a  higher  likelihood  of  belonging  to 
the  line.  Hence  the  method  is  much  more  robust  than  searching  for  fp’s  of  the  grout  line 
intersections. 

The  robustness  of  the  hough  line  method  is  observed  in  the  following  where  a  rectangular  paper 
is  imaged  with  it’s  corners  removed.  Hough  line  transform  has  no  issue  with  determining  the 
virtual  corners  of  the  rectangle  that  would  otherwise  have  been  present. 


Figure  3-32  Hough  line  transform  and  infinite  line  intersection  applied  to  a  rectangle  with 
corners  missing. 
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4.  SGL  Experimental  Apparatus 

4.1  Overview 

As  described  in  the  previous  chapter,  simulation  of  the  SGL  with  the  SA  is  limited  due  to  the 
difficultly  of  justifying  the  propagation  model  in  a  multipath  environment.  Essentially  the 
performance  results  of  the  SA  are  readily  predictable  given  a  specific  multipath  model. 
However,  to  access  the  actual  performance  of  the  implemented  SA  based  SGL  it  is  necessary 
to  have  a  set  of  statistical  measurements. 

Measurements  based  on  spread  spectrum  modulated  signals  will  be  made  with  chipping  rates 
of  up  to  20  Mcps.  From  these  the  channel  impulse  response  (CIR)  can  be  measured.  Having 
up  to  four  transmitter  antennas  in  various  locations,  the  simultaneous  multipath  propagation 
channels  can  be  represented  which  are  difficult  to  model  justifiably.  The  other  main  objective  of 
the  measurement  system  is  that  of  the  LOS  detection  and  LOS  processing  which  get’s  to  the 
main  contribution  of  the  SA. 

A  wireless  system  consisting  of  several  transmit  antennas  and  a  receiver  that  can  be  moved 
while  taking  measurements  has  been  created  and  is  used  for  the  experimental  component  of 
this  contract.  This  wireless  system  can  emulate  the  transmission,  reception  and  processing  of 
the  LTE  downlink  signal  reference  signals  that  the  SGL.  The  carrier  frequency  can  be  tuned 
within  the  operating  range  of  700  MHz  to  2.2  GHz  with  various  modulation  options  from  a  simple 
tone  to  a  full  20  MHz  bandwidth  modulation.  In  phase  III  the  objective  will  be  to  expand  this 
bandwidth  further  to  100  MHz  that  will  account  for  the  evolution  of  the  LTE  standard. 
Modifications  can  be  made  to  the  existing  hardware  system  to  provide  this.  This  will  involve  the 
receiver  and  transmitter  using  a  step  frequency  such  that  the  100  MHz  is  covered  by  multiple 
bands  each  of  20  MHz  instantaneous  bandwidth.  Figure  5.1  shows  a  typical  configuration  in  the 
experimental  room  of  ENA401  at  the  University  of  Calgary.  The  idea  is  to  locate  the  4  antennas 
such  that  they  have  a  variable  amount  of  multipath  and  in  some  cases,  as  the  trajectory  is 
changed  that  the  propagation  will  vary  from  LOS  to  NLOS. 
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Figure  4-1  Floor  plan  layout  of  the  LTE  propagation  room  ENA401  at  the  University  of  Calgary 
The  wireless  measurement  system  that  was  developed  will  be  described  in  section  4.2.  Section 
4.3  describes  the  processing  of  the  single  receiver  channel  with  a  reference  antenna  channel. 
Section  4.4  describes  the  multi-antenna  setup  and  the  overall  test  environment.  Section  4.5 
provides  experimental  results  and  Section  4.6  an  overall  summary  and  some  conclusions. 

4.2  Wireless  measurement  system 

The  wireless  measurement  system  developed  essentially  consists  of  a  set  of  four  AN’s  and  a 
UE.  Each  AN  consists  of  an  RF  amplifier  and  antenna  that  is  mounted  to  the  wall  in  the  test 
area.  The  UE  consists  of  a  small  antenna  and  LNA  for  receiving  the  AN  transmitted  signal.  It  is 
also  bundled  with  a  webcam  for  the  associated  CV  measurements.  The  transmit  signal  is 
generated  by  a  National  Instrument  software  radio  transceiver,  denoted  as  the  USRP.  A  similar 
matching  USRP  is  used  in  the  receiver  to  facilitate  the  processing.  A  block  diagram  of  the 
overall  architecture  of  the  system  is  given  in  Figure  4.2. 
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Figure  4.2  Overall  block  diagram  of  the  LTESA  hardware  system 

The  USRP’s  for  the  transmitter  and  receiver  functions  are  programmed  via  LabView.  This 
provides  maximum  flexibility  for  the  various  experimental  setups  that  must  be  implemented.  As 
LabView  is  a  rather  cumbersome  software  development  environment  with  discovered 
limitations,  it  will  only  be  used  for  the  signal  modulation  for  the  transmitter  and  the  sample 
extraction  in  the  receiver.  The  SA  processing  as  well  as  the  general  SGL  processing  will  be 
done  off  line  via  Matlab  and  C.  In  Figure  4.1  it  is  observed  that  the  receiver  has  two  antennas. 
One  is  used  as  an  optional  reference  antenna  and  the  other  antenna  is  main  SA  receiver 
antenna.  The  receiver  has  a  separate  laptop  that  is  used  to  run  labview  to  collect  the  data  from 
the  USRP  receiver  and  to  collect  the  CV  data.  The  transmitter  is  run  from  a  separate  laptop  and 
USRP.  The  reason  for  this  is  that  the  receive  and  the  transmitter  can  be  physically  separated 
without  the  inconvenience  of  having  to  tie  them  together  with  data  connections.  Also  there  is  a 
10  MHz  synchronizing  signal  that  is  used  to  frequency  lock  the  receiver  and  the  transmitter. 
This  is  so  that  the  frequency  instability  of  the  receiver  can  be  removed  in  some  experiments  to 
study  impacts.  An  additional  laptop  may  be  required  for  the  CV  processing  of  the  webcam 
output.  The  issue  is  that  if  a  single  PC  services  boththe  webcam  and  the  USRP  then  the 
connections  are  a  bottleneck  resulting  in  issues  of  irregular  sampling  such  that  the  camera  and 
wireless  data  are  not  easily  synchronized.  Additional  components  for  the  system  are  the 
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precision  turntable  which  is  used  to  move  the  antenna  as  an  alternative  to  hand  motion.  The 
precision  turntable  motion  implies  that  the  webcam  is  not  needed  as  the  antenna  trajectory  is 
easily  calibrated.  This  turntable  is  driven  by  the  PC  used  for  the  receiver.  Additionally,  the 
receiver  PC  controls  the  selection  of  up  to  four  transmit  antennas  as  shown  in  Figure  4.3. 
Shown  is  the  Digital  Input  Output  (DIO)  interface  block  connected  to  the  receiver  PC  via  a  USB 
that  provides  the  output  digital  signals  for  sequencing  the  antenna,  turntable  stepper  motor  and 
receives  inputs  from  the  trigger  switch. 


Figure  4-3  Block  diagram  of  the  antenna  selection,  trigger  and  turntable 
A  block  diagram  of  the  USRP  is  given  in  Figure  4.4.  It  consists  of  a  transmitter  and  a  receiver 
channel  that  can  be  switched  to  the  antenna  ports.  In  the  transmitter,  generic  quadrature 
modulation  is  used  with  a  complex  sampling  rate  of  up  to  33  MHz.  Likewise  in  the  receiver 
quadrature  demodulation  is  used  with  a  sampling  rate  of  up  to  33  MHz  and  16  bits  of 
quantization.  A  gigabit  Ethernet  connection  is  used  for  all  of  the  data  exchange  between  the 
host  PC  and  the  USRP.  This  is  a  limitation  as  the  full  modulation  bandwidth  at  16  bits  I  and  Q  in 
addition  to  the  Ethernet  protocol  overhead  exceeds  the  1  Gbps  capacity  of  the  connection. 
Hence  if  higher  IQ  sampling  rates  are  used  then  it  is  necessary  to  limit  the  sampling  block  size. 
The  USRP  has  internal  buffering  such  that  high  rate  burst  sampling  can  be  accommodated 
provided  that  the  average  data  traffic  required  is  within  the  1  Gbps  limitation.  The  limitations  of 
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the  Ethernet  connectivity  is  also  the  reason  why  two  PC’s  are  used,  one  dedicated  for  the 
transmitter  and  the  other  for  the  receiver. 


Figure  4.4  Block  diagram  of  the  USRP 

The  Labview  programming  of  the  receiver  USRP  is  a  simple  state  machine.  Unfortunately  due 
to  the  buffering  of  the  data  it  is  not  possible  to  precisely  control  measurement  segments  of  the 
USRP.  Rather  it  samples  for  the  4.5  seconds  at  a  rate  of  200,000  samples  per  second  which 
are  then  transferred  to  the  hard-drive  of  the  receiver  PC.  Hence  the  turntable  and  receiver 
antenna  sampling  are  driven  by  a  separate  process.  The  USRP  and  turntable/antenna  switch 
are  run  as  two  separate  processing  threads.  The  states  of  the  LabView  program  for  gathering 
the  USRP  data  is  simply  an  initialization  phase,  sampling  phase,  data  storage  phase  and  a 
termination  phase  that  are  run  in  sequence.  An  issue  with  the  USRP  is  that  every  new 
measurement  interval  implies  a  new  initiation  with  causes  the  LO  to  re-lock.  This  uncontrolled 
phase  cannot  be  tolerated  in  the  SA. 

Prior  to  starting  the  sampling,  the  frequency  needs  to  be  adjusted.  If  the  rx  is  offset  by  lOOppm 

then  the  result  is  1575.42AfflzxlO  5  =16 kHz .  This  offset  is  reduced  by  adding  an  offset  to  the 
receiver  signal.  It  is  reasonable  to  assume  that  this  can  be  hardwired  in  bringing  the  uncertainty 
down  to  about  1  kHz.  The  minimum  sampling  rate  is  about  200  kHz.  Hence  there  is  obviously 
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no  issue  with  aliasing  and  the  input  signal  can  be  low  pass  filtered  to  4  kHz  and  decimated  by  a 
factor  of  20  to  10  kHz  effective  sampling  rate  prior  to  storage. 

There  are  various  hard  wired  variables  used  as  listed  in  the  following  table: 


variable 

description 

samp_rate 

Sampling  rate  set  to  the  minimum  of  200k  for  a  CW  signal  and  up  to  33 

MHz  for  a  spread  spectrum  signal 

fLO 

LO  frequency  set  to  1.57542  GHz  for  the  GPS  related  sampling  but  can 

be  anywhere  from  700  MHz  to  2.2  GHz 

channel 

RX2  is  used  primarily. 

Number  of 

samples 

Can  be  set  arbitrarily  high.  Typical  for  the  CW  SA  analysis  a  million 

samples  are  used  for  a  5  second  SA  epoch. 

A  brief  description  of  the  LabView  programs  will  be  given.  The  LV  panel  for  the  stepper 
for  the  turntable  as  well  as  the  receiver  antenna  is  given  in  Figure  5.5.  Note  that  it  is 
necessary  to  run  these  operations  from  a  separated  processing  thread  so  that  it  does 
not  interfere  with  the  USRP  sampling  operations.  Settable  controls  are  the  number  of 
events  which  are  essentially  steps  of  the  turntable  stepper.  Every  step  is  set  at  8  msec 
and  after  four  steps  (32  msec)  the  antenna  is  switched.  The  turntable  direction  can  be 
toggled  as  clockwise  or  counter  clockwise.  Finally  the  turntable  can  be  disabled  such 
that  only  antenna  toggling  is  used.  Figure  5.6  gives  the  LV  block  diagram  which 
consists  of  a  single  for  loop  that  is  incremented  every  8  msec.  A  formula  block  is  used 
to  give  the  digital  outputs  to  the  DIO  for  the  antenna  switch  and  the  stepper  sequencing. 


Figure  4.5  Panel  of  stepper/RF  receiver  antenna  switch  control 
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Figure  4.6  LV  block  diagram  of  the  stepper/antenna  switch  control 
The  formula  for  the  stepper  sequencing  is  given  by 
X2*(mod(X1 ,4)+4*mod(X1 ,4)+4+8*X3) 

Here  XI  is  the  integer  index  of  the  loop,  X2  is  the  boolean  variable  enabling  the  turntable  to 
move  and  X3  is  the  direction.  The  four  least  significant  bits  are  used  by  the  stepper  driver. 

The  panel  for  the  USRP  control  is  shown  in  Figure  5.7.  There  are  no  user  controls  to  set.  The 
panel  only  indicates  when  the  sampling  is  taking  place  and  the  coerced  USRP  variables.  Also 
shown  is  the  display  of  the  samples  which  provides  an  indication  of  the  quality  of  the  sampling. 
In  this  instance,  the  sampling  of  the  reference  and  the  receiver  antenna  can  be  seen  in  the 
graph.  The  reference  should  be  the  first  to  start  and  should  be  fairly  uniform  in  amplitude.  The 
receiver  antenna  in  this  case  shows  a  multipath  faded  signal  as  the  antenna  is  moved. 
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Figure  4.7  USRP  Control  panel 
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Figure  4.8  USRP  Control  block  diagram  for  the  initialization  state. 

Figure  4.8  shows  part  of  the  initialization  of  the  USRP  with  the  hardwired  variables  of 
LO  frequency,  receiver  gain  in  dB  and  sampling  rate  as  well  as  the  receiver  channel 
RX2.  The  actual  sampling  occurs  in  this  state.  Figure  4.9  shows  the  state  for  writing 
the  I  and  Q  waveforms  to  files  on  the  hard  drive. 


h  stop  4 


Figure  4.9  USRP  Control  block  diagram  for  the  initialization 


Figure  4-10  shows  the  block  diagram  of  the  two  receiver  antennas,  switch  and  turn  table 
stepper.  Figure  4-1 1  is  a  picture  of  the  receiver  antenna  on  the  turntable. 
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Figure  4-10  Block  diagram  of  antenna  switch  and  stepper 
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Figure  4-11  Photo  of  receiver  antenna  and  reference  antenna  setup  on  turntable 


Figure  4-12  shows  a  close-up  of  the  receiver  antenna  and  the  web-cam.  The  antenna 
is  magnetically  mounted  so  that  it  can  be  moved  such  that  it’s  phase  center  can 
coincide  with  the  webcam  image  center. 
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Figure  4-12  Close-up  of  the  receiver  antenna  and  the  web-cam 


4.3  Processing  of  received  wireless  signal 

As  discussed  in  the  previous  section  the  USRP  samples  the  I  and  Q  signals  of  the  antenna 
output  as  a  contiguous  block.  The  reference  antenna  and  the  receiver  antenna  outputs  are 
commutated  by  a  switch  into  the  USRP.  The  file  of  typically  800k  samples  will  be  stored  on  disk 
and  then  read  into  Matlab  for  offline  processing.  This  processing  will  be  described  in  this 
section.  The  complete  Matlab  program  is  listed  at  the  end  of  the  chapter. 

Initially  the  IQ  data  is  read  into  the  matlab  workspace  and  combined  into  a  set  of  complex 
samples  r.  A  plot  of  |r|is  shown  in  Figure  4-13.  The  transitions  are  due  to  the  switching 
between  the  two  antennas  which  in  the  present  configuration  occurs  every  32  msec.  The  slow 
variation  is  due  to  the  antenna  being  moved  by  the  turntable  in  an  environment  where  there  is 
multipath. 
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Figure  4-13  Plot  of  |  r  |  showing  the  variations  of  the  two  multiplexed  antenna  responses 
The  first  task  is  to  determine  the  decomposition  of  the  samples  of  r  such  that  the  reference  and 
received  channels  can  be  separated.  While  the  timing  is  not  perfect,  it  is  known  that  the  dwell 
time  of  the  sampling  is  32  msec  for  each  antenna.  Typically,  as  in  the  example  shown  in  5-1 3, 

|  r  |  will  be  quite  different  for  the  two  antennas.  Hence  the  processing  is  to  slide  a  template 
across  with  intervals  of  32  msec  of  samples  each.  The  best  alignment  is  when  the  variance  of 
|  r  |  is  the  smallest  in  each  interval.  Note  that  if  the  intervals  contain  the  switch  transition  then 
|  r  |  will  change  abruptly  within  the  interval  resulting  in  a  large  variance.  The  average  interval 
variance  is  denoted  by  G  and  is  plotted  in  Figure  4-14  as  a  function  of  the  offset  of  the  template 
pattern.  There  is  a  well  defined  minimum  which  is  then  used  to  give  the  offset  required. 
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Figure  4-14  Plot  of  the  average  interval  variance  of  |  r  |  as  a  function  of  the  offset 
With  this  offset  the  data  can  then  be  segregated  into  the  two  antenna  responses.  In  the  current 
example  which  is  plotted  in  Figure  4-15  the  red  samples  are  for  the  receiver  antenna  and  the 
green  samples  are  for  the  reference  antenna.  Note  that  the  green  reference  samples  show 
less  variation  as  expected. 


Figure  4-15  Plot  of  segregated  data 

As  it  is  undesirable  to  have  to  control  the  timing  of  the  sampling  it  may  not  be  obvious  which  is 
the  reference  channel.  However,  this  is  generally  clear  from  the  variations  of  |  r  |  Also  for 
whatever  reason  (to  be  investigated)  the  first  interval  is  sometimes  too  short  and  hence  there  is 
a  switching  transition  that  occurs  within  the  interval.  This  results  in  an  interval  of  samples  of 
high  variance  that  should  be  discarded.  This  is  shown  in  Figure  5-16.  Figure  5-17  shows  a  plot 
of  the  real  and  imag  components  of  each  sample  after  they  have  been  derotated  in  phase 
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based  on  an  approximate  estimate  of  the  offset  frequency.  Note  that  only  the  reference  signal 
samples  are  used  for  the  estimate  of  the  offset  frequency. 
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Figure  4-16  Sample  of  incomplete  interval  having  a  high  variance  of  |  r  \ 


Figure  4-17  Typical  example  of  a  FFT  of  the  sample  segments  indicating  the  estimation  of  the 
offset  frequency 
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Figure  4-1 8  Plot  of  the  real  and  imaginary  components  of  each  sample.  Green  and  blue 
correspond  to  the  reference  frequency  and  red  and  black  to  the  receiver  components 
Initially  an  approximate  phase  shift  between  the  reference  and  the  receiver  channels  can  be 
determined  by  multiplying  the  signal  in  a  reference  channel  by  the  conjugate  of  the  interval  of 
samples  in  the  adjacent  receiver  channel.  The  angle  of  the  complex  valued  signal  provides  an 
estimate  of  the  difference  phase  that  unfolds  as  the  antenna  is  moved.  Unwrapping  is  applied  to 
the  phase  angle  of  these  samples.  The  phase  trajectory  is  plotted  in  Figure  5-19. 
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Figure  4-19  Approximate  phase  excursion  as  the  antenna  is  moved. 
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Clearly  the  phase  excursion  in  Figure  5-19  is  not  overly  accurate  but  it  does  provide  the 
experimenter  a  quick  assessment  of  the  quality  of  the  data.  Other  checks  can  be  done  such 
as  stability  of  the  signal  s  when  the  turntable  is  not  moved.  This  is  shown  in  Figure  4-20. 


time  (sec) 

Figure  4-20  Response  when  the  antenna  is  held  constant  position 

Fig  4-21  shows  the  partitioning  of  the  sampling  with  the  red  and  green  as  the  portions  of  the 
samples  that  for  the  reference  and  receiver  channels.  The  blue  samples  are  discarded  as  they 
correspond  to  the  transition  regions. 
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Fig  4-21  Detail  of  Figure  4-20 


4.4  Simultaneous  CV  and  LTE  outputs  and  processing 
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Unfortunately,  the  USRP  operation  consumes  significant  bandwidth  of  the  USB  interface  to  the 
laptop.  Likewise  the  webcam  demands  a  near  maximum  transfer  rate  when  the  webcam  is 
streaming  video  at  the  rate  of  more  than  20  frames  per  second.  The  Windows  operating  system 
cannot  accommodate  two  simultaneous  USB  streaming  sources  and  therefore  if  the  units  are 
run  at  the  same  time  then  the  operation  will  be  jerky  with  unpredictable  delays.  While  the  CV 
can  accommodate  interruptions  to  handle  the  USRP  operation,  the  USRP  needs  to  be  reset 
every  time  which  causes  excessive  delays  and  the  phase  coherency  will  not  be  maintained. 

Therefore  it  is  necessary  to  use  two  laptops  to  accommodate  both  the  webcam  input  and  the 
USRP.  The  labtops  are  time  synchronized  with  a  common  HW  switch  that  is  triggered  at  the 
start  and  end  of  the  sampling.  Once  the  raw  samples  have  been  stored  they  will  be  processed 
in  parallel  on  the  two  laptops.  An  Ethernet  connection  between  the  two  computers  combines 
the  processed  CV  and  USRP  data  on  a  single  computer.  This  is  then  fed  into  the  overall  SLAM 
algorithm. 

At  this  stage,  due  to  the  interfacing  issue,  it  is  not  possible  to  run  the  HSA  hardware  completely 
in  real  time.  However,  the  processing  of  the  raw  samples  is  relatively  fast  such  that  the  SLAM 
processing  starts  essentially  as  soon  as  the  sampling  is  completed. 
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Fig  4-22  Overall  flow  of  the  data  sampling  and  processing  algorithm 
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4.5  Upgraded  Experimental  system  based  on  Nl  RF  PXI 

An  upgraded  indoor  localization  experimental  system  based  on  wireless  RF  signal 
measurements  is  established.  The  hardware  of  the  system  mainly  consists  of  Nl  RF  PXI 
system,  Nl  USRP,  antennae,  turntable,  and  switches.  The  software,  based  on  LABVIEW, 
implements  the  configuration  of  the  Nl  PXI  system  and  USRP,  control  of  the  external  devices, 
and  synchronization  of  the  overall  system.  RF  signal  measurements  are  acquired  and  streamed 
to  the  PXI  system,  where  they  are  stored  and  can  be  used  to  apply  the  data  post-processing. 
The  following  resources  are  needed  for  the  system  setup: 

■  Hardware 


-  Nl  PXIe-1075  (18-Slot  3U  PXI  Express  Chassis  with  AC  -  Up  to  4  GB/s) 

-  Nl  PXI-8 135  (PXI  Express  Embedded  Controller) 

-  Nl  PXIe-5663E  (6.6  GHz  Vector  Signal  Analyzer  (VSA)  With  RF  List  Mode) 


-  Nl  PXIe-6555  (200  MHz  Digital  Waveform  Generator,  High  Speed  Digital  10  (HSDIO)) 
with  Nl  CB/SCB  -  2162  single-ended  digital  I/O  accessory 

Generally,  we  call  the  combination  of  the  above  devices  as  PXI  system,  which  is  shown  in 
figure  4-23. 


Figure  4-23:  PXI  system 


-  Nl  USRP-2920  (50  MHz  to  2.2  GHz  Software  Radio) 


-  Turntable  and  controller 
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-  Dipole  LTE  Antenna 

-  Coaxial  SP4T  switch 
■  Software 

-  Lab  VIEW  2012 

-  MathWorks  MATLAB  2013a 


The  application  of  this  experimental  system  is  to  do  indoor  wireless  localization  based  on 
synthetic  array  (SA).  Basically,  in  the  experiment  we  need  a  receiver  moving  along  a  known 
trajectory  while  receiving  wireless  signals  from  access  nodes  (AN)  and  find  its  position  relative 
to  the  AN. 
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Figure  4-24:  Overview  of  the  indoor  localization  experimental  system 

Figure  24  shows  the  overview  of  the  experimental  system.  The  experimental  room  is  at  ENA401 
at  the  University  of  Calgary.  Basically,  a  transmitter  is  placed  in  a  turntable  moving  in  a  circle.  4 
receiving  antennae  are  located  at  each  corner  of  the  room  to  receive  signal  from  the  transmitter. 
Then  the  4  received  signals  will  be  scanned  by  a  switch  alternatively  and  fed  to  a  receiver, 
where  the  signal  is  proceed  to  obtain  the  measurements  for  localization. 


125/143 


Figure  24  also  shows  the  implementation  of  the  system.  Specifically,  the  receiver  is 
implemented  by  USRP.  Figure  3  shows  the  USRP  and  the  transmitting  antenna  on 


Figure  4-25:  URSP  and  transmitting  antenna  on  the  turntable 


the  turntable.  Shown  in  figure  26  is  the  receiving  antenna  in  the  corner  of  the  room. 
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Figure  4-26:  Receiving  antenna  placed  in  the  corner  of  the  room 


Figure  4-27:  Switch  used  to  scan  signals  from  receiving  antennae 

The  switch  is  controlled  by  HSDIO  and  the  accessory  board  as  shown  in  figure  4-28. 
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Figure  4-28:  HSDIO  and  accessory  board  to  control  switch 

The  receiver  is  implemented  using  VSA  in  PXI  system,  as  shown  in  figure  4-23. 
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Figure  4-29:  Signal  processing  functional  blocks  of  the  experimental  system 
Figure  4-29  describes  the  signal  processing  functional  blocks  of  the  experimental  system. 
Generally,  a  bandlimited  pseudorandom-noise  (PN)  code  is  generated  and  stored  in  a  file  using 
MATLAB.  USRP  reads  the  file  to  fetch  baseband  waveform,  and  up-converts,  interpolates,  low 
pass  filters,  and  mixes  the  signal  to  RF  for  transmitting.  In  the  receiver  part,  the  VSA  receives 
the  RF  signal,  then  down-converts,  band  pass  filters  and  digitizes  the  signal.  Decimation  is  then 
applied  to  the  signal  and  the  received  baseband  waveform  is  stored  in  a  file  for  further 
processing.  In  the  post-processing,  the  received  signal  is  dispread  using  the  local  PN  code, 
from  which  the  amplitude,  phase,  time  of  arrival  (TOA)  measurements  can  be  obtained  for  the 
localization  algorithm. 
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Unlike  the  devices  realised  by  Nl  which  can  be  control  by  LABVIEW  directly,  the  way  to  control 
the  switch  and  turntable  need  to  be  given  more  thoughts,  so  as  to  make  the  system  integration 
easier  and  more  efficient.  The  HSDIO  provides  flexible  software  define  way  to  generate  digital 
waveform.  Along  with  accessory  board,  HSDIO  can  control  the  switch  through  LABVIEW.  Figure 
8  shows  the  HSDIO  accessory  board  and  switch  as  well  as  the  connection 


Figure  4-30:  HSDIO  accessory  board  and  switch 

between  them.  For  the  turntable,  its  controller  can  be  programed  using  LABVIEW  based  driver, 


Figure  4-31:  Turntable  LABVIEW  based  drive 
as  shown  in  figure  4-31 . 


To  sum  up,  all  the  hardware  involved  in  the  system  can  be  control  by  LABVIEW,  which  makes 
the  integration  of  the  system  using  one  single  host  program  possible.  Figure  4-32  shows  the 
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Figure  4-32:  Overview  of  the  block  diagram  of  the  host  program 

overview  of  the  block  diagram  of  the  host  program.  Figure  1 1  shows  the  Front  panel  of  the  host 
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Figure  4-33:  Front  panel  of  the  host  program 

program.  In  the  experiment,  we  just  need  to  configure  the  setting  of  each  module  in  the  program 
front  panel,  and  simply  press  start  button.  Then  the  turntable,  switch  and  receiver  will  start 
working  in  a  synchronized  way,  and  the  received  baseband  waveform  will  be  recorded  in  a  file 
specified  by  the  program  for  further  processing. 

FREQUENCY  SYNCHRONIZATION  OF  TRANSMITTER,  RECEIVER  AND  SWITCH 
CONTROLLER 

Frequency  synchronization  plays  an  essential  role  in  the  integration  of  such  a  system  involving 
multiple  hardwires.  Figure  4-34  explain  the  frequency  synchronization  scheme  of  the  system. 
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Figure  4-34:  Frequency  synchronization  scheme  of  the  system 

As  shown,  USRP,  VSA  and  HDSIO  share  a  same  reference  clock  of  10MHz  in  VSA,  and 
generate  the  respective  system  clock  using  respective  PLL.  Then  each  system  clocks  will  be 
used  to  generate  other  clocks  to  drive  the  modules  in  this  device.  The  synchronization  of 
transmitter  and  receiver  in  a  hardwired  connection  way  is  an  option  to  make  the  synchronization 
easier  in  the  first  place,  and  will  be  replaced  by  some  other  algorithm  in  the  receiver  later. 


To  implement  above  scheme,  features  of  VSA,  USRP  and  HSDIO  are  utilized  in  the  design. 
Generally,  VSA  can  export  internal  reference  clock  to  the  ports  in  its  panel.  As  shown  in  figure 
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Figure  4-35:  VSA  reference  clock  outputs  and  configuration 

13,  through  the  configuration  in  the  program,  the  internal  reference  clock  of  VSA  are  available  in 
the  2  ports  highlighted  in  the  figure.  USRP  and  HSDIO  can  accept  external  clock  source  as 
reference  clock,  as  shown  in  figure  14  and  figure  15.  So  we  can  simply  connect  their  reference 


Figure  4-36:  USRP  reference  clock  input  and  configuration 
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Figure  4-37:  HSDIO  reference  clock  input  and  configuration 

clock  input  port  to  the  clock  source  port  in  VSA  panel  and  finish  the  setting  for  frequency 
synchronization. 

EVENT  SYNCHRONIZATION  OF  RECEIVER  AND  SWITCH  CONTROLLER 
Event  synchronization  of  the  receiver  and  switch  controller  must  be  performed  to  identify  the 
signal  from  different  channels.  Figure  38  gives  the  event  synchronization  scheme. 


Figure  16:  Event  synchronization  scheme  of  receiver  and  switch  controller 
Basically,  HSDIO  uses  its  sample  clock  to  generate  synchronized  switch  control  signal  and 
event  trigger  signal.  While  the  switch  control  signal  is  fed  to  the  switch,  the  event  trigger  signal 
is  routed  to  the  VSA  to  trigger  VSA  start  recording  data. 


133/143 


Features  of  HSDIO  and  VSA  are  used  to  implement  the  scheme  above.  HSDIO  is  able  to 
generate  digital  I/O  waveforms  and  events  signal  (called  marker  event)  in  specified  patterns  by 
script,  as  shown  in  figure  39. 


Script 

script  MyScript  a 

repeat  forever 

Generate  switch_ctrl_wav  markerO  (0) 
end  repeat 
end  script 


Figure  4-39:  Configuration  of  HSDIO  to  generate  digital  I/O  waveform  and  event  trigger 
signal 

The  digital  I/O  waveforms  then  will  be  exported  to  I/O  at  the  accessory  board  to  control  the 
switch,  while  event  signal  will  be  export  to  the  port  in  HSDIO  panel  shown  in  figure  4-40. 
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Figure  4-40:  HSDIO  event  signal  output  and  configuration 

Since  VSA  is  able  to  accept  external  event  trigger  input,  as  shown  in  figure  4-41 ,  to  trigger  the 
data  recording,  we  can  simply  connect  this  input  port  to  the  event  trigger  source  in  HDSIO  port 
and  finish  the  event  synchronization  setting. 


Figure  4-41:  VSA  event  trigger  input  and  configuration 

Figure  20  shows  the  generated  switch  control  signal  and  VSA  event  trigger  signal  capture  by 
logic  analyzer.  As  can  be  seen,  the  VSA  event  trigger  signal  can  mark  the  beginning  of  the 
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Figure  4-42:  Switch  control  signal  and  VSA  event  trigger  signal  generated 
channel  0  data  recording.  Figure  21  shows  the  baseband  waveforms  recorded  in  file.  Obviously, 
we  can  separate  the  date  of  each  channel  very  easily  because  the  beginning  of  the  date  record 
corresponds  to  that  of  the  first  channel. 
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Figure  4-42:  Waveform  record  of  4  channels 


Experimental  results 

Figure  4-43  shows  the  measurements  and  location  results  of  an  experiment.  In  this  experiment, 
only  1  receiving  antenna  is  used  and  is  placed  in  the  position  highlighted  by  the  red  star,  while 
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Figure  4-43:  Waveform  record  of  4  channels 

the  transmitter  is  moved  along  a  circle  represented  by  green  circle  in  right  graph.  The  phase 
measurements,  after  removing  circle  clip,  are  shown  in  the  left  graph  with  green  dots.  We  then 
apply  localization  algorithm  to  phase  measurements,  and  the  blue  circles  in  the  right  graph 
show  how  the  position  estimate  converges,  which  tells  that  the  position  estimated  by  receiver  is 
very  close  to  the  actual  position. 
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5.  Discussion 


The  initial  objective  of  this  research  was  to  develop  elements  of  a  system  that  could  use  the  LTE 
wireless  signals  of  opportunity  for  self  geo-location  (SGL)  of  a  handheld  navigation  device 
(HND).  This  implied  that  the  HND  was  not  registered  with  the  wireless  system.  The  issues  of 
not  being  able  to  assume  that  the  LTE  sources  were  in  known  locations  or  time  synchronized 
implied  that  some  form  of  SLAM  processing  was  necessary.  It  was  determined  that  provided 
that  it  could  be  assumed  that  the  LTE  access  nodes  were  stationary  and  frequency  synchronized 
(not  code  or  offset  synchronized)  then  SLAM  processing  would  be  able  to  simultaneously 
determine  where  the  LTE  access  nodes  are  and  their  offsets  in  time  delay  as  well  as  the 
trajectory  of  the  HND  relative  to  the  LTE  network  of  access  points.  SLAM  also  provided  the 
most  general  framework  for  implementing  the  synthetic  array  component.  Extensive  details  of 
the  SLAM  algorithm  implementation  and  analysis  is  contained  in  part  B  of  this  final  report. 

While  the  multipath  of  the  LTE  emissions  can  be  modelled  such  that  the  SLAM  algorithm  can  to 
some  degree  average  out  the  variances  of  the  errors  of  the  HND  trajectory,  this  is  shown  to  be 
the  difficult  component  of  achieving  a  sub-meter  accuracy.  Hence  CV  observables  were  added 
into  the  mix  of  measurement  constraints.  These  were  shown  to  be  significant  in  terms  of  the 
overall  SLAM  in  that  CV  positioning  is  highly  accurate  over  short  trajectory  lengths  but  subject 
to  drift.  Contrary  wireless  with  multipath  is  rather  inaccurate  due  to  multipath  but  not  subject 
to  the  same  drift  as  with  the  CV.  Hence  CV  and  wireless  observables  are  complementary.  Of 
course  other  observables  can  be  used  such  as  inertial  but  this  was  beyond  the  scope.  Inertial 
will  have  the  same  issues  as  CV  in  terms  of  drift.  However,  as  inertial  drifts  with  time  and 
spatial  movements,  CV  only  drifts  as  a  result  of  nonlperfect  calibration  of  the  camera,  image 
noise,  non-ideal  features  and  imager  discretization.  CV  is  not  affected  by  time  lapse.  Again 
inertial  can  be  construed  as  being  complementary  to  both  CV  and  wireless.  It  should  therefore 
be  included  in  the  mix  of  observables. 

An  emerging  development  that  we  did  not  have  bandwidth  to  cover  is  that  of  cooperative 
positioning  (CP)  which  fits  very  nicely  within  the  overall  SLAM  framework.  This  is  the  joint 
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location  of  several  mobile  nodes  simultaneously  which  share  information.  CP  requires  some 
form  of  messaging  between  the  mobile  nodes  or  messaging  to  a  central  node.  It  is  easily 
demonstrated  in  the  SLAM  formation  given  in  part  B  that  CP  can  significantly  multiply  the 
positioning  information  that  is  mutually  available  reducing  the  variance  of  the  desired  variables 
associated  with  the  trajectory  of  the  mobile  units. 

Future  work  would  focus  on  the  combining  of  the  three  sensor  types,  wireless,  CV  and  inertial 
into  a  single  instrumented  experimental  testbed.  While  experimental  apparatus  was  developed 
for  the  wireless  and  CV,  these  need  to  be  combined  in  a  more  unified  way  such  that  the  relative 
significance  of  the  measurement  modalities  can  be  assessed. 
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